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1. INTRODUCTION 



In a recent report [1] a mixed numerical/analytical approach to the computation of a ring- 
symmetric spacecraft exhaust plume was presented. The numerical scheme had been implemented in 
a code named "JET" which is capable of generating whole-plume flow fields, while the analytic 
approximation is restricted to the ring- symmetric centered rarefaction waves (CRW) that flank the 
plume. The present report is intended to serve as a supplement to (1) in providing details on the 
computational scheme and the code JET. 

The spacecraft exhaust flow (Fig. 1 of (11 ) is idealized as a ring- symmetric steady isentropic 
expansion of an ideal gas. The nozzle lips are assumed sharp; the supersonic flow from the exit 
surface of the ring-nozzle is assumed uniform, and the background is considered to be perfect 
vacuum. 

The standard scheme for computing such idealized ring-plumes is the classical (direct) method of 
characteristics [2] . At a preliminary phase of the present laser exhaust study, a code AXSYM [3] was 
WTitten for computing ring-plumes using this method. A notorious shortcoming of the direct method 
of characteristics is that the solution grid is highly irregular, being formed by the (oblique) 
intersection of the C and the C ~ families of characteristic lines. We first encountered a dilTiculty 
with this grid while seeking a scheme for integrating the molecular opacity along a straight line [1] . 
This computation would have required rather complex coding for the geometiy of intersection 
between a straight line and an irregular grid. It seemed preferable to opt for a computation scheme 
that would produce a more regular grid, even at the e.xpense of some loss of accuracy. Such scheme 
is the inverse marching method of characteristics [4] . 

Generally the marching in this type of scheme is in the downstream direction, i.e., the y direction 
in our case. Grid points are located on a succession of constant-y rows, thereby introducing a 
measure of regularity in the solution grid. Just two rows have to be stored in the computer core 
memory - the "old" row and the "new" row, whereas in the direct method of characteristics whole 
grid-image matrices are required to reside simultaneously in core memory. 

The first version of the JET code was based on the inverse marching scheme given by Zucrow and 
Hoffman (Section 12-5 in [4J ), where the flow variables were the two cartesian velocity components. 
The computation seemed accurate everv-where, except within the centered rarefaction wave (CRW). 
In an attempt to replicate a planar CRW (Prandtl- Meyer flow), the numerical solution exhibited an 
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instability : Mach number increased along the (low pressure) boundary characteristic line, rather than 
remain constant. 

A qualitative explanation for this instability is the following. Flow gradients in a CRW are 
inversely proportional to distance from the comer, so that the inverse marching scheme gives rise to 
an amplification of interpolation errors at every marching step, leading to an apparently divergent 
(unstable) numerical solution. Increasing the order of interpolation from linear to cubic did not 
eliminate the instability. 

Looking for a scheme that would replicate a planar CRW accurately, we tried the modified 
marching idea as presented by Zucrow and Hoffman for 1-D time-dependent flows (Sections 19-6(a) 
and I9-6(j) of [4] ). In this scheme new grid points are determined by forwardly extending a "primary" 
family of continuous characteristic lines from old grid points. The primary family in a CRW is the 
characteristics fanning out from the corner (we assume it is the C"^ family). By choosing this 
modified scheme, the interpolation for trace points obtained from reversely extended C lines was 
eliminated. However, the corresponding interpolation for the transverse C~ characteristics 
remained, and with it the aforementioned instability. 

In order to replicate a planar CRW, we had to replace the flow variables by the Riemann 
invariants (v±6). In a C"*" planar CRW, the Riemann invariant (v + 0) is uniformly constant, so 
that the interpolation in (v + 0) due to reversely extending C ~ characteristics introduces no error at 
all. This scheme, which we named SIMA (Semi Inverse Marching Algorithm), was indeed verified to 
replicate a planar CRW exactly, when implemented in the code JET. 

The plan of this report is the following. In Ch. 2 we supplement the description of the numerical 
scheme given in Ch. 2 of [I] , by adding more details on the computational procedure. A description 
of the code JET is given in Ch. 3, and the code listing is reproduced in Ch. 4. 

Note on symmetry’ : 

The code JET has two symmetry options. When DELTA = I a ring-symmetry is in effect; when 
DELT.“\ = 0, a planar symmetiy’ is in effect. .A.n axisymmetric jet exiting in the y direction from the 
same nozzle aperture along the X axis can readily be computed by replacing all terms in the code that 
correspond to sin(0)/y in the compatibility equations (2.1-1), by cos(0)/x. In that case the coding is 
virtually unchanged, and the only care that should be exercised is for the difference equations for new 
grid points on or near the y axis. Also, all reference to the analytic approximation of the ring- 
symmetric CRW [1] should be deleted in this case, as it is designed specifically for ring-symmetry. 



2. THE COMPUTATIONAL SCHEME 



A basic description of the senii inverse (SIMA) and inverse marching schemes was given in Ch. 2 
of [1] . We supplement this description by specifying the slightly modified definition of Riemann 
invariants in the code, and by giving information about some ancillary computations. 



2.1 Riemann Invariants 

The compatibility equations whose integration constitutes the numerical solution to the governing 
equations [1] are expressed in terms of the Riemann invariants as follows : 

Along C (V - 0)4 = (V - 9)2 + CO sin^24 sin024 An / >^24 

( 2 . 1 - 1 ) 

Along C“ (v + 9)^ = (v + 9)j + (0 sinpj^ sin0j4 A^ / y,4 

The Riemann invariants (v±9) are modified for convenience, by adding a constant to both v and 
0. The new definitions of v(M) and 9 are : 

v(M) = — r arctan(Fq)+ arctan(q) 

q=(M"-l)’‘^^ (2.1-2) 

9 0 - 0L 

Thus, in a Prandtl-Meyer flow with entry .Vlach number of M|, the modified values of both v(M) 
and 0 vanish as M-»go. As a consequence, in a C"*” Prandtl-Meyer flow the modified invariant 
(v + 9) vanishes uniformly. In this modified form, the computation of M from v(M) is readily done 
by performing standard Newton-Raphson iterations (in RFUNC), using the derivative : 

v'(q) = -(r -l)[(l + r q*)(l + q*)] (2.1-3) 
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2.2 The Integration Scheme for a New Grid Point 



The integration scheme has been sketched in Ch. 2 of [1] . It is performed in INVMAR for inverse 
marching points or in SEMINV for semi-inverse marching (SIMA) points. The computational 
scheme is specified via the following seven-step procedure : 



INVMAR (Inverse Marching) 

(a) Grid : At this stage the new grid point has already been defined. 

(b) Predictor : Flow variables are the interpolated (linear nearest-neighbor) value on the old row 
for a point having the new grid X coordinate (x_^). 

(c) Centered variables : Denote the Riemann invariants by 



RM = (v + 0) 
RP = (v-0) 



( 2 . 2 - 1 ) 



then centered values for segments (1,4) and (2,4) (using code notation) are : 
RM14 = (RMl + RM4)/2 RP14 = (RPl + RP4)/2 

RM24 = (RM2 + RM4)/2 RP24 = (RP2 + RP4)/2 



( 2 . 2 - 2 ) 



All other centered flow variables are computed from the centered Riemann invariants by calling 
RFUNC. 

(d) Inverse Extension : old trace points Xj, x, are evaluated from the geometrical relations 

Along C “ .... = (x^ - x,) tan(0„ - (i,,) 

(2.2-3) 

Along C + .... y„^„ - y,„ = (x., - x^) tan(0jj + 

(e) Interpolation : find Riemann invariants RM, RP at old trace points Xj and X 2 through nearest- 
neighbor linear interpolation by calling INTERP. 

(f) Integration : Using the compatibility relations in finite-difierence form (2.1-1) with segment- 
centered coefficients, compute iteration-updated values of Riemann invariants at new grid point. 
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(g) Corrector : if values of Riemann invariants and old trace points Xj, X 2 are not sufficiently 
convergent, resume the procedure at step (c) above. 

SEMINV (Semi Inverse Marching - SIMA) 

(a) Grid : New grid point (x^) is determined as part of the SIMA scheme at step (d) below. 

(b) Predictor : Flow variables are those of point 

(c) Centered variables : Identical to step (c) above. 

(d) Semi-Inverse Extension : new grid point x^ and old trace point Xj are evaluated from the 
geometrical relations in Eq. (2.2-3) above. 

(e) Interpolation : find Riemann invariants RM, RP at old trace point X| through nearest- 
neighbor linear interpolation by calling INTERP. 

(f) Integration : Identical to step (f) above. 

(g) Corrector : Identical to step (g) above, except for replacing X 2 in the convergence test by x_^. 

2.3 Boundary Conditions 

On the vacuum side the boundarv^ conditions (p = 0) can only be approximately implemented in a 
method of characteristics scheme. We do so by ending the computation on a certain "final" C fan 
characteristic line that starts out with a sufficiently high Mach number Mj. at the corner (typically 
Mp= 34). The marching computation of new grid points on the boundary' C characteristic via the 
SIMA scheme is identical to that of C"*" characteristics within the ring-symmetric CRW. It is noted 
that under this boundary scheme some outflow takes place through the boundaiy characteristic line, 
so that the total mass flow through a row y= y^^^^ decreases slightly as increases. 

At the nozzle exit the boundary conditions are assumed to be uniform outflow in the radial (y) 
direction with Mach number M,. At the nozzle lip. the SIMA integration starts out from a presumed 
planar CRW (Prandtl-Meyer flow) at the corner (i.e., the associate CRW in the terminology of Ch. 3 
in [1] ). 

At the plane of symmetry (x=0) the boundary condition is simply 0 = rt/2. However, this 
condition is implemented indirectly, by assuming that the flow at virtual grid points with x < 0 is a 
mirror-image of the flow at the corresponding X > 0 points. The reason is that when a new grid 
point of Xj=0 or of Xj sufficiently close to zero is considered for inverse-marching integration, the 
inversely extended trace point (Xj.y^ju) can be at X < 0. Considering the subtraction of 6^^ from 0 as 
in Eq.(2.1-2), the reflection rules are ; 
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RM RP + (71-20^) 



(2.3-1) 



RP -* RM - (71-20^) 

where values on the left and right of the -» symbol correspond to values left and right of x=0. This 
boundary condition is implemented in INTERP. 



2.4 Continuum Breakdown Surface 



As an informative option, the code JET can compute (in PLUMES) points on a surface of 
continuum breakdown [5,6,7] , which is defined as a line of constant B, where B is given by : 



B= -(u/(p)p-‘(dp/dS) 
- 1/2 

<p = 4(7iy) <r n a 



(2.4-1) 



When the standard isentropic relations for p and n in terms of M are substituted in (2.4-1), the 
fiow speed is expressed as U = Ma and the streamwise gradient of M is expressed in cartesian 
coordinates, we get : 

B = X^(7tY/8)'^^M‘ [i+((Y-1)/2)mY^^'^**^ ' [m^cos0 +M^sin0] 

(2.4-2) 

^0 = (2 n ^) 

Note that the sign of B has been chosen as positive for expansion flows. This definition is 
preferred to taking an absolute value of the flow gradient, since it assures proper interpolation of B 
even if its spatial distribution goes through B = 0. 



Due to the dependence of B on a spatial gradient, its numerical evaluation (see BREAK) is 
attributed to mid-grid points both in X and in y. 
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3. THE JET CODE 



In this chapter we provide a concise description of the JET code according to its version at the 
time of the JET018 run. This description is intended as an aid in reading the code listing which is 
given in Ch. 4. 

The plan of this chapter is as follows. Array variables that constitute the mainstay of the 
computational scheme are described in Section 3.1. Auxiliary array variables that are used primarily 
for processing the information generated by the numerical scheme, are described in Section 3.2, 
followed in Section 3.3 by a list of major parameters that control the computation (some of them also 
serve as run data). Finally, all subroutines are listed and described in Section 3.4. 



3.1 Main Variables 

The array variables used for the computational scheme are organized in two labeled CO.VIMON 
groups. The first group /VECS/ is designed to hold two grid rows - the old row designated by suffix 
F and the new row designated by suffix N. The second group /CHAR.-\C/ are characteristic-indexed 
arrays that hold information about continuous characteristic lines. This characteristic information is 
used in two ways : it is incorporated in the SIMA computational scheme for the CRVV region, and it 
is used to store data for optional plotting of characteristic lines (see PLUMES and PRINT). 

The basic organization is that the new arrays (suffix N) are those in which values are stored during 
the course of the marching computational procedure. At the end of each marching step, values are 
transferred from new arrays to old arrays (suffix F); this is done in MOVE. In the array listing 
below, we indicate in parenthesis the subroutine (or subroutines) in which that new array is defined. 

/VECS/ • 



XN(I) 


X coordinate of grid point 1. 


(GRIDN) 




RMN(I) 


modified Riemann invariant 


(v + 0) at grid point I. 


(BEGIN. INVMAR. 




LOADC). 






RPN(l) 


modified Riemann invariant 


(v — 0) at grid point I. 


(BEGIN. INVMAR. 




LOADC). 






.MN(I) 


Mach number at grid point I 


(BEGIN, INVMAR, LOADC). 


MUN(I) 


Mach angle p at grid point I. 


(BEGIN, INVMAR, LOADC). 


TETAN(I) 


true (unmodified) flow angle 9 at grid point I. (BEGIN, INVMAR, LOADC). 
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BN(I) 



value of breakdown parameter B at point 1-1/2 (and at half a marching step 
back in y as well). (BREAK). 

XTEMP(I) used for auxiliary computation of 1-1/2 grid points in PLUMES. 



/CHARAC/ 

XCHARN(KC) 

YCHARN(KC) 

RMCARN(KC) 

RPCARN(KC) 

TCHARN(KC) 

MUCARN'(KC) 

CSIGNN'(KC) 

MCHARN(KC) 

MCHARl(KC) 



X coordinate of point on characteristic line number KC. (BEGIN, 
SEMINV, PLUMES). 

y coordinate of point on characteristic line number KC. (BEGIN, SEMINV, 
PLUMES). 

modified Riemann invariant (v + 0) of point on characteristic line number 
KC. (BEGIN, SEMINV). 

modified Riemann invariant (v - 0) of point on characteristic line number 
KC. (BEGIN, SEMINV). 

true (unmodified) flow angle 0 at point on characteristic line number KC. 
(BEGIN, SEMINV). 

Mach angle ft at point on characteristic line number KC. (BEGIN, 

SEMINV). 

sign of characteristic line number KC. It has value 1 for C"^ and value - 1 
for C ~ . Note that upon reflection of a C ^ line from the symmetry plane 
(x= 0), the sign value is changed from 1 to - 1. (BEGIN, SEMINV). 

Mach number at point on characteristic line number KC. (BEGIN, 

SEMINV). 

Mach number at Prandtl-Meyer's fan characteristic number KC at the 
corner. It is defined initially and is not changed during the run. (BEGIN). 



3.2 Auxiliary Variables 

In addition to the major arrays mentioned above, there are several groups of auxiliary arrays that 
do not affect the computational scheme, but are intended for informative processing of the results. 
These groups are /PLUME/, /IPLUME/, /THICKY/, /THICKX/, /GRP/. /PLUME/ is used to 
preserve points on special lines for later plotting (in a separate code). /THICKY/ and /THICKX/ 
are for storing values of radial (y) and lateral (x) molecular opacities. The group /GRP/ is used in 
conjunction with comparative computation of the ring-symmetric CRW flow according to the analytic 
approximation [1] . 



8 



/PLUME/ (PLUMES, PRINT) 



XPL(J.IPL) 

YPL(J.IPL) 


X coordinate at marching step J of special line number IPL. 
y coordinate at marching step J of special line number IPL. 



/IPLUME/ (PLUMES, PRINT) 

KPL number of special lines computed in PLUMES. 

ITYPL(IPL) index indicating the type of special line number IPL. 

/THICKY/ (OPACY, PRINT) 

XTH(J) X coordinate on boundary characteristic line at marching step J, from which 



TH(J) 


radial opacity is integrated. 

radial opacity computed by y-integration from the boundary point defined by 
XTH(J) (up to current YN). 



/THICKX/ (OPACX, PLUMES, PRINT) 

YXI(JXl) y coordinate of printed row number JXI (the index JXI counts just rows that 



Xl(UXl) 


have been printed). The row to be printed next upon calling PRINT is the row 
having YF near YXI(JXl). 

lateral (x) molecular opacity [1] at point XF(I), for printed row JXI. It is 
obtained by numerically integrating the solution obtained from the JET 
computation (see OPACX). 


XlP.M(LJXl) 


same as XI(I.JXI) e.xcept that the Prandtl- .Meyer solution is used to estimate 
the flow at grid points XF(I). 


X1GRP(I,JX1) 


same as above, except that the analytic approximation to a ring-symmetric 
CRW [1] is used to estimate the How at grid points XF(I). 


X1APP(1,JX1) 


same as XIGRP(I,JXI) except that the numerical integration is replaced by an 
approximate closed-form expression [IJ . 


X1F(1,JX1) 


stores grid points XF(I) of printed row JXI. 



/GRP/ (PRINT, HMSET, MFUNC, HINTER, MATCH) 

DMINV increment of inverse Mach number for array .VIHINV(I). 



.MHINV(I) 


inverse Mach number array (from 0 to l/.MEXIT), from which the H(M) 
function can be evaluated (HVISET). 


HMV(I) 


values of the H(M) function evaluated by numerical integration. It is used to 
compute this function by interpolation. (HMSET, HINTER). 
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3.3 Major Parameters 



Parameters that define and control a particular run (such as the maximum y for the marching 
computation, the number of grid points on a row and many more) are defined in INIDAT. (The 
code JET has no input file and no READ statements). The major control parameters are grouped in 
/PAR/ (floating point) and in /IPAR/ (integers); thermodynamic data are grouped in /STAG/. 

We indicate in the listing the subroutines in which the labeled COMMON group or a particular 
parameter is defined (or sometimes referred to). 



/PAR/ (INIDAT) 

MEXIT nozzle exit Mach number (Mj). 

MEIN Mach number of the final (boundary) CRW characteristic at the comer (Mj.). 

YMAX maximum value of y for the marching scheme. When YF.GE.YMAX the run is 
terminated. 



DYO 

DY 

DYNEXT 

STAB 

DELTA 

PSIl 

PSIF 

SIGxMA 

FR.ACG 

EPSIL 

TETLI.M 

TETSYM 



initial marching step, 
current marching step, 
next marching step (YSTEP). 

stability coefficient for marching step (STAB.LE.l). (See YSTEP). 
symmetry index. DELTA = 0 for plane symmetry; DELTA= 1 for ring-symmetry, 
angle of Prandtl-Meyer fan characteristic at exit conditions (measured from x 
axis). 

angle of final (boundary) Prandtl-Meyer fan characteristic, 
collision cross-section (<T). 

the number of intervals initially allocated to the CRW fan is a FRACG fraction of 

the total number of intervals (KFO-1). (see BEGIN). 

convergence parameter (small number). (INVMAR, SEMINV). 

flow angle (from x axis) of the limiting (p= 0) velocity vector of the flow at the lip- 

centered Prandtl-Meyer fan. 

PA1-2*TETLI.M for reflection transformation (see INTERP). 



/IPAR/ (INIDAT) 

JMAX maximum number of marching steps. If J.GE.JMAX run is terminated. 

KFO initial (and maximum) number of grid points in a row. 

KF current number of grid points in the old row. 

KN current number of grid points in the new row. 
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ITERO 


maximum number of iterations for the integration of the compatibility relations 
(see INVMAR and SEMINV; also used in RFUNC, PLUMES). 


IM, IP 
J 


search indices for interpolation subroutine INTERP. (see INVMAR, SEMINV). 
current row index (also index of a marching step). 


KF2 


defined as 2*KF; not used in present version. 



IDEL, JDEL increments for printing grid point I and row J (see PRINT). 



JYXI 


number of rows to be printed in a run. 


JXI 


index of printing row, to be printed next (see PRINT). 


I LEAD 


index I at the first grid point on current new row, where the SIMA integration 
commences. Initially this point corresponds to the leading characteristic of the 
CRW. (see GRIDN, BEGIN). 


ILEADF 

KCLEAD 


value of I LEAD for current old row. 

index in the characteristic array for the characteristic line that corresponds to the 
new grid point 1 = ILEAD (see GRIDN). Initially KCLEAD = 1. 



/STAG/ (INIDAT) 

RHOO, NO stagnation density and number density. 

PO, TO, AO stagnation pressure, temperature and sound speed. 

MDOTl mass flow rate from ring-nozzle (only from the x > 0 halO. (See PRINT). 
/ICHARA/ (BEGIN) 



KCHARP 


number of C"^ characteristic lines for which data is stored (either for SIMA 


KCHARM 


computation or for subsequent plotting). 

number of C~ characteristic lines for which data is stored (only for subsequent 
plotting). 


KCHARO 


total number of characteristics for which data is stored, i.e., 
KCHAR0= KCHARP + KCHARM. 



II 



3.4 Description of JET subroutines 



MAIN PROGRAM 

The main program performs two functions. The first section (up to statement 1) is the initial set 
up; it is performed just once. The second section is the marching loop with the step index J. This 
program can be read as a flow chart of the overall computational procedure. 

INIDAT is for setting up run data. In BEGIN the initial conditions for the marching 
computation are set up. A single marching step is performed by calling MARCH, and the loading of 
new row vectors into old row vectors is done by calling MOVE. The call to YSTEP is for the first 
computed marching step. All remaining calls are for informative tasks (see HMSET, BREAK, 
OPACY, PLUMES, PRINT). Run is terminated when either YF.GE.YMAX or when J.GE.JMAX. 

NOTE ON EXEC : The only special feature in the EXEC is retaining the output unit 7 file for 
optional post-plotting. The printed output (unit 6) is the system's standard (default). 



INIDAT 

Initial data definition and preliminary data computations. The data is defined by statements rather 
than by reading an input file. The meaning of major parameters was described in Section 3.3 above. 
User is invited to modify the data definitions, particularly of run-control parameters such as YMAX, 
JYXl and YXI(JXI) (for printing JYXI selected rows). 



BEGIN 

Here all initial values (prior to beginning of marching schemes integration) are loaded into all 
major computational arrays (Section 3.1). Also, values of the key integer parameters KCHARP, 
KCHARM, KCHARO, I LEAD, KCLEAD and KF are defined. 

In the first loop (loop 1) we define an initial family of C"^ characteristic lines for the lip-centered 
CRW, by storing the Mach number of the Prandtl-Meyer fan characteristics in the array 
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MCHARI(KC). Note that the fan characteristics are generated at equal RP intervals, since the flow 
variables are RM and RP. However, a different division might also be acceptable. 

The next step is the deflnition of initial values for all characteristic arrays, first the C arrays 
(loop 2), then the C arrays (loop 21). The C~ characteristic lines are needed just for informative 
output (post-plotting), so the present version contains just one C ~ line. The user may modify that. 

The remaining grid points (altogether KFO grid points are initially available) are uniformly 
distributed across the nozzle opening, and the row arrays are loaded \vith the corresponding nozzle- 
exit flow variables (loop 3). 



PRINT 

The main task of this subroutine is the printing of flow variables at grid points of selected rows. 
The printing of a row is selected when YF is close to a predefined array YXI(JXI). Following the 
printing, JXI is updated by adding 1. 

For comparison, additional flow variables are printed for each row. These are computed from the 
analytic approximation to a ring-symmetric CRW [1] , by calling MATCH. Also, lateral molecular 
opacities of various kinds of approximation are computed by calling OPACX, and are printed for 
each grid point within the CRW. 

Following the row printing (statement 120), arrays intended for post-processing (plotting of special 
lines) are printed and subsequently written on output unit 7. This is done once per run, just before 
run termination. 



FIN 

This subroutine is called when an error is encountered, in order to terminate the run. Note that 
the run is terminated by deliberately introducing an error of computing SQRT(-l), which is done in 
order to trigger the printing of calling sequences by the operating system. 
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MARCH 



This subroutine performs a single marching step by calling the proper computational subroutines 
at an appropriate sequence. It can be read as a flow chart of the entire computational scheme. First 
the segment of the new row suitable for SIMA computation is calculated by calling SEMINV. Then 
new grid points for that part of the new row for which inverse marching integration is to be 
performed, are generated by calling GRIDN. The results of the SEMINV computation, which were 
stored in characteristic arrays, are now loaded into row arrays by calling LOADC. Finally, the 
computation of the new row is completed by calling INVMAR which computes the flow at the 
remaining grid points by the inverse marching scheme. 



INVMAR 

This is one of the two central subroutines for computing the flow at new grid points (the other is 
SEMINV). Here the inverse marching scheme is used. The computational procedure follows the 
seven- step description given in Section 2.2 above. Note that the initial value of the search indices IM 
and IP is not redefined at each call to INTERP, since it is assumed that IM and IP do not change 
much at consecutive calls to INTERP, so that search efficiency is enhanced by not starting the search 
from an arbitrary point (such as either end of the row). 



SEMINV 

This is the subroutine performing the SIMA scheme for computing the flow at new grid points 
located along continuous characteristic lines of the lip-centered CRW (at prescribed y-marching 
steps). The essence of the computational procedure of this subroutine was given as a seven-step 
description in Section 2.2. The same remark about IM given in the preceding INVMAR description 
applies here as well. 

The main loop ( 100) is over all characteristic lines, including some C ~ lines in addition to the C ^ 
lines. Thus, the array CSIGNF(KC) is used to get the appropriate expressions for either C"*" or C~ 
characteristics. It is noted that while normally the characteristic segments through points 1 and 2 are 
C and C respectively, this is reversed when a C ~ rather than a C ^ line is computed via the 
SIMA scheme. In this case, which is characterized by having CSIGNF(KC).LT.O, the Riemann 
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invariants integrated along segments (1,4) and (2,4) are interchanged. This is done in the few 
statements just preceding and following statement 21. 

An additional capability of this subroutine is to treat a change of a C characteristic line into a 
C” line upon reflection from the symmetry plane (x=0). This is done by first computing a new grid 
point having X4.LT.0, and then changing its sign after setting CSIGNN(KC)= — 1 (statements just 
preceding statement 30). It is also possible to skip the computation of a particular characteristic by 
setting CSIGNN(KC) = 0. This feature is not exploited in the present version. 

Finally, we note that not all characteristic lines computed here are part of the marching flow 
computation. Only those with indices KC between KCLEAD and KCHARP are. All other 
characteristic lines are computed just for informative purposes (post-plotting). 



RFUNC 

Here M, MU, TETA are computed from the two Riemann invariants RM, RP. The computation 
of M is performed by a Newton- Raphson iteration using Equations (2.1-2) and (2.1-3) given in 
Section 2.1 above. 



INTERP 

This subroutine starts by fmding through a search procedure the grid interval (I, I+l) that 
contains a given point X. Then the Riemann invariants are computed for this point by linear 
interpolation, and returned in RM, RP. Note that X may be negative, which accounts for the 
relatively elaborate search logic in the determination of 1, and for the reflection transformation (as in 
Eq. (2.3-1) above) preceding the last two statements of the subroutine. 



INTERX 

This interpolation routine performs an inverse task to that of INTERP. in that it finds the point 
XO that corresponds to a given linearly interpolated value of the flow variable VARO. It is used in 
PLUMES to compute the location of a breakdown surface point on a new row of x-centered and y- 
centered grid points 
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BREAK 



This subroutine computes the new breakdown parameter array BN(I). The computation is based 
on the description given in Section 2.4 above. 



OPACY 

Here the radial (Y) molecular opacity array TH(J) is computed. At each marching step J, a new 
boundary grid point XTH(J) is added, then the radial opacities at all preceding boundary points are 
updated by adding the contribution of the gas layer between the current old and new rows. Note that 
since grid points on adjacent rows are not located on equal-x columns, this procedure requires X- 
interpolation by calling INTERP. 



PLUMES 

This is a user-defined subroutine, where up to 10 special lines can be computed and subsequently 
retained on output unit 7 for post-processing (plotting). The type of the line ITYPL(IPL) and a 
parameter VPL(IPL) that defines the line, are computed through user-inserted statements in the 
section preceding statement 2000. Then an additional point on the current new row is computed for 
each line type. The available types are clearly stated in comments. Note that characteristic lines 
have already been computed in SEMINV using the SIMA scheme, regardless of whether they are part 
of the solution grid to the flow field, or are just computed for informative purpose. It is the user's 
choice which of these lines (if any) are to be saved in the /PLUME/ arrays for subsequent post- 
processing (plotting). 



GRIDN 

This subroutine computes the grid points in that segment of the new row for which the flow is 
computed by the inverse marching scheme (in INVMAR). Initially, this segment extends from x=0 
to the new row grid point which lies on the leading characteristic of the lip-centered CRW. However, 
since the leading characteristic is reflected from the symmetry plane (x=0) at some point, this 
segment steadily shrinks in size as the marching proceeds. The remedy is to declare the next-to-the- 
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leading characteristic line (KC=2) as the beginning of the segment for SIMA integration, by setting 
KCLEAD = 2. This process of increasing KCLEAD is repeated whenever it is deemed necessary. 
The criterion in the present version for the minimal KCLEAD is that the inverse-marching segment 
should be at least twice DXl - the average CRW grid interval (loop 1, the two statements following 
DXl = ...). Also, ILEAD is redefined for each row according to XLEAD/DXl + 2 in order to achieve 
a row of relatively uniform grid intervals throughout. The result is that the number of grid points in 
a row is initially KFO, but eventually it decreases due to both increase of KCLEAD and decrease of 
ILEAD. 

YSTEP 

In this subroutine the next marching step DYNEXT is computed at the end of the current 
marching step. It is defined as the smallest step obtained by foru'ard intersection of C ~ and C 
characteristics from adjacent grid points. Note that the actual value of DYNEXT is reduced by a 
"stability" factor STAB, and that DY is also limited by the growth-rate factor DDY and by DYMAX 
(see MAIN PROGRAM). 



MOVE 

Here old row arrays (loop 1) and old characteristic arrays (loop 2) are loaded with values of flow 
variables from corresponding new arrays, in preparation for the next marching step. As a result of 
this organizational feature, informative computations (e.g. BREAK, OPACY) that require both new 
and old rows, have to be performed prior to calling MOVE. 



OPACX 

Here lateral (X) opacities that correspond to the number of expected collisions of a fast molecule 
invading the CRW in the —X direction, are computed. All opacities, except XIAPP(I), are computed 
by numerical integration. In loop 1 we compute the opacity contribution of the segment lying just 
outside the computational boundary characteristic (MFIN), assuming a Prandtl- Meyer flow. This 
additional opacity is denoted XIO. If the flow is ring-symmetric, XIO is recalculated using the analytic 
approximation [1] to estimate the flow field at the fringes of the ring-synunetric CRW (see also the 
closed form expression for t in [1] ). 
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The computation of opacity arrays starts after statement 14. First, the opacity at each grid point 
is set to XIO. Thus, even though the numerical flow computation does not include the fluid outside 
the boundary characteristic line, the opacity integration includes an estimate of that "missing" part, 
i.e., of XIO. In typical case computations of a ring-symmetric CRW [1] we found that the maximum 
value of XIO w'as about 0.16., which indicated that as far as interaction with invading ambient 
molecules is concerned, the approximation MFIN= 34 was a reasonable substitute for MFIN= oo. 

The next step is the computation by numerical integration of three approximations to the lateral 
opacity : XI(I,JXI), XIPM(I,JX1), XIGRP(I,JX1). (Note that when the flow is ring-symmetric, the 
approximation XIPM(I,JXI) obtained by assuming a Prandtl-Meyer flow is usually grossly 
exaggerated). The opacity XIGRP(I,JXI) is based on the analytic approximation to a ring-symmetric 
CRW [1] , and is reasonably close to XI(I,JXI) which is obtained from the numerical solution to the 
flow field. Finally, a simplified closed-form integration of lateral opacity [1] is computed as 
XIAPP(I,JXI) (loop 3). Thus, the quantitative difference between XI(I,JXI) and XIGRP(I,JXI) is an 
indication to the degree of accuracy achieved by the analytic approximation to a ring-symmetric 
CRW [1] , while the difference between XIGRP(I,JXI) and XIAPP(I,JXI) indicates the level of error 
introduced by the closed-form integration of lateral opacity [1] . 



LOADC 

Here flow variables of new grid points computed via the SIMA scheme (SEMINV) are loaded into 
new row' arrays from corresponding characteristic arrays. 



NUFUNC 

This function computes the modified v(M) value as given by Eq. (2.1-2). Note that presently 
NU0=0(see INTDAT). 



HMSET 

This subroutine is called just once from the .MAIN PROGRAM. Its task is to set up the arrays in 
/GRP/, so that the function H(M) [1] can be evaluated by interpolation (in HINTER). There is also 
an informative printout of various derivatives (see Ch. 3 of [1] ) generated in this subroutine. 
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MFUNC 



This subroutine is called by HMSET in order to compute functions of Mach number that serve in 
the computation of H(M). The output variable F is the integrand for the integration leading to 
H(M). 

HINTER 

This subroutine computes H(M) by linear interpolation in inverse Mach number, using the /GRP/ 
arrays computed in HMSET. 



MATCH 

This subroutine is called from PRINT to compute the Mach number according to the analytic 
approximation of a ring-symmetric CRW [1] , for point (YF,XF(I)). MOB is the associate Mach 
number M(0,P), which is preserved in the array MCHARI(KC) for all CRW characteristics that are 
used in the SIMA computation. Hence the Vlach number M(a,P), denoted by MAB can be 
computed directly from the analytic approximation [1] to the area function at (YF,XF(I)) by calling 
AREAF. Since typically M(O.P) is not knowm, we also compute the Mach number via the inverse- 
problem procedure (11 , denoting the resulting Mach numbers by suffix I : MOBI for M(O.p) and 
MABI for M(a,P). The inverse-problem iterative procedure [Ij is performed in loop 1. resulting in 
MOBI. From MOBI the value of MABI is computed through the area function approximation as for 
MAB above. 



AREAF 

This subroutine computes the Mach number M that corresponds to the area function F (Eq. 
(3.2-1) of [1] ). The computation is done by Newton- Raphson iterations, and it has been found to 
converge when M.GT.l (and when M — 1 is not much smaller than 1). 
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4. THE JET CODE LISTING 



C#0PTI0NS LIST 
C JET018 

C "JET" A SEMI-INVERSE MARCHING CHARACTERISTICS METHOD FOR RING JETS. 
C USING RIEMANN INVARIANTS RM=(NU+TETA) , RP=(NU-TETA) AS FIELD 
C VARIABLES. 

IMPLICIT REALX8(A-H,L-Z,$) 

REALX4 XPL,YPL 

COMMON /PLUME/XPL(1002,10),YPL(I002) 

COMMON /IPLUME/KPL,ITYPL(10) 

COMMON /VECS/XF(I01),RMF(I0I),RPFa01),MF(101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3 TETAN(101),BN(101),XTEMP(I0I) 

COMMON/THICKY/XTH(l 002), TH( 1002) 

REALXA YXI,XI,XIPM,XIGRP,XIAPP,XIF 

COMMON /THICKX/YXI(20),XI(101,20),XIPM(101,20),XIGRP(101,20) 

I ,XIAPP(I01,20),XIF(I01,20) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,GI0,G1I,G12,GI3,G14,GI5 
1 G16,G17,GI8,G19,G20 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RH00,N0,P0,T0,A0,MD0T1 
COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

1 KF2,IDEL,JDEL, JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON /ROW/YF,YN,DXF,DXN 

COMMON /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), 

1 RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), 

2 TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 

4 MCHARK92) 

COMMON /ICHARA/KCHARP, KCHARM, KCHARO 
COMMON /GRP/DMINV,MHINV(101),HMV(101) 

COMMON /IGRP/KHM 
C 

PRINT 101 
101 FORMAT('l') 

J = 1 

IF(J.EQ.l) STOP 
CALL INIDAT 
PRINT 101 
CALL HMSET 
PRINT 101 
CALL BEGIN 
CALL MARCH 
CALL OPACY 
CALL PLUMES 
CALL PRINT 
J=2 

CALL PLUMES 
CALL MOVE 
CALL OPACY 
CALL PRINT 
CALL YSTEP 
1 J=J+1 

C DY WAS DETERMINED BY THE PREVIOUS CALL TO GRIDN. 

DY=DMIN1( DYNEXT, DY)«DDY, DYMAX) 

C INTEGRATE BY ONE Y-STEP 
CALL MARCH 

C BREAKDOWN PARAMETER (BF(D). 

CALL BREAK 

C SPECIALLY DESIGNATED LINES (FOR PLOTTING). 

CALL PLUMES 

C STORE NEW LINE (N) IN OLD LINE (F). 

CALL MOVE 

C COMPUTE RADIAL MOLECULAR OPACITIES 
CALL OPACY 

C Y-STEP IS VARIABLE, SO JMAX IS USED AS END-OF-RUN CRITERION. 

IF(YF.GE.YMAX) JMAX=J 
C PRINT FIELD AT MOST RECENT Y. 

CALL PRINT 
C NEXT Y-STEP. 



JETOOOl 

JET0002 

JET0003 

JET0004 

JET0005 

JET0006 

JET0007 

JET0008 

JET0009 

JETOOlO 

JETOOll 

JET0012 

JET0013 

JET0014 

JET0015 

JET0016 

JET0017 

,JET0018 

JET0019 

JET0020 

JET0021 

JET0022 

JET0023 

JET0024 

JET0025 

JET0026 

JET0027 

JET0028 

JET0029 

JET0030 

JET0031 

JET0032 

JET0033 

JET0034 

JET0035 

JET0036 

JET0037 

JET0038 

JET0039 

JET0040 

JET0041 

JET0042 

JET0043 

JET0044 

JET0045 

JET0046 

JET0047 

JET0048 

JET0049 

JET0050 

JET0051 

JET0052 

JET0053 

JET0054 

JET0055 

JET0056 

JET0057 

JET0058 

JET0059 

JET0060 

JET0061 

JET0062 

JET0063 

JET0064 

JET0065 

JET0066 

JET0067 

JET0068 

JET0069 

JET0070 

JET0071 

JET0072 
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. JETPR 



FORTRAN A1 



CALL YSTEP 
IF(J.LT.JMAX) GO TO 
STOP 
END 






JET0073 
JET007A 
JET0075 
JET 131176 



C 

C 



SUBROUTINE INIDAT 
SUBROUTINE NUMBER 1 

IMPLICIT REALX8(A-H,L-Z,$) 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3 TETAN(101),BN(101),XTEMP(101) 

■ COMMON/THICKY/XTH( 1002) ,TH( 1002) 

REALX4 YXI,XI,XIPM,XIGRP,XIAPP,XIF 

COMMON /THICKX/YXI(20),XI(101,20),XIPM(101,20),XIGRP(101,20) 

1 ,XIAPP(101,20),XIF(101,20) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14, 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI , PAI2, DEG, XC, YC, MEXIT, MFIN, YMAX, DYO , DY, DYNEXT , 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 
COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

1 KF2,IDEL,JDEL, JYXI, JXI,ILEAD,ILEADF,KCLEAD 

COMMON /ROW/YF,YN,DXF,DXN 

PAI=4.D0*DATAN(1 .DO) 

PAI2=2.D0XDATAN(1.D0) 

DEG=180.D0/PAI 

AR=8.3143D3 

AV=6.022D 26 

AW=7 .27D0 

RHO0=0.0075D0 

T0=2300.D0 

G=1.54D0 

D=2.5D-10 

MEXIT=4.D0 

MFIN=34.D0 

XC=0.5D0 

YC=2.5D0 

DELTA=0 CORRESPONDS TO PLANE SYMMETRY 
DELTA=1 CORRESPONDS TO CYLINDRICAL SYMMETRY 
DELTA=1.D0 
FRACG=0 .6D0 
EPSIL=l.D-8 
ITER0=20 
KF0=101 
JMAX=1001 
STAB=0.50D0 
DDY=1 . 05D0 
DYMAX=0.5D0 
YMAX=50.D0 
DY0=YC/250.D0 
IDEL=1 
JDEL=1 

POINTS FOR PRINTING FLOW FIELD AT YF=YXI(JXI) 

JXI = 1 
JYXI=11 
DYXI=5.D0 
YXI(1)=YC+0.5D0 
YXI(2)=YXI(1)+2.D0 
10 = 2 

DO 1 1=10, JYXI 

YXI(I)=YXI(IO)+DYXI*DFLOAT(I-IO) 

CONTINUE 

IF(KFO.GT.lOl) CALL FIN(lOl) 

IF(JMAX.GT.lOOl) CALL FIN(102) 

IFCFRACG.GT.l.DO .OR. FRACG.LT.O.) CALL FIN(103) 

IF( JYXI .GT.20) CALL FIN(104) 

IF(DELTAX(1. DO-DELTA). NE.O.) CALL FINCIOS) 

NO=RHOOXAV/AM 

A0=DSQRT(GXARXT0/AM) 

PO=AR*RHOOXTO/AW 



JET0077 
JET0078 
JET0079 
JET0080 
JET0081 
JET0082 
JET0083 
JET0084 
JET0085 
JET0086 
JET0087 
G15, JET0088 
JET0089 
JET0090 
JET0091 
JET0092 
JET0093 
JET0094 
JET0095 
JET0096 
JET0097 
JET0098 
JET0099 
JETOlOO 
JETOlOl 
JET0102 
JET0103 
JET0104 
JET0105 
JET0106 
JET0107 
JET0108 
JET0109 
JETOllO 
JETOlll 
JET0112 
JET0113 
JET0114 
JET0115 
JET0116 
JET0117 
JET0118 
JET0119 
JET0120 
JET0121 
JET0122 
JET0123 
JET0124 
JET0125 
JET0126 
JET0127 
JET0128 
JET0129 
JET0130 
JET0131 
JET0132 
JET0133 
JET0134 
JET0135 
JET0136 
JET0137 
JET0138 
JET0139 
JET0140 
JET0141 
JET0142 
JET0143 
JET0144 
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SIGMA=PAIXDXX2 

LAMDA0=1.D0/(DSQRT(2.D0)XSIGMAXN0) 

Gl=(G-l.DO)/2.DO 

G2=(G+1 .DO)/(2.DOX(G-l.DO)) 

G3=G/2.D0 

GA=(G+l.DO)/(G-l.DO) 

G5=DSQRT((G+1.D0)/(G-1.D0)) 

G6=1.D0/(G-1.D0) 

G7=2.D0/(G+1.D0) 

G8 = (0.5D0*(G+1.D0)xx2/(G-1.D0))*3((1.D0/(G+1.D0))* 

1 ((G+l.DO)/(G-l .DO))*X((G-l.DO)/(G+l.DO)) 

G9=(G+3.D0)/(2.D0X(G-1.D0)) 

G10=(7.D0-3.D0XG)/(2.D0X(G-1 .DO)) 
G11=(2.D0/(G+1.D0))XX(1.D0/(G-1.D0)) 
G12=DSQRT((G+1.D0)/(G-1.D0))-1.D0 
G13=(2.D0-G)/(2.D0X(G-1.D0)) 

G1A=G/(2.D0*(G-1.D0)) 

G15=(G+1.D0)/(3.D0-G) 

G16=(G+1.D0)/A.D0 

G20=LAMDA0XDSQRT(PAIXG/8.D0) 

ZETAl=G5XDATAN(DSQRT(MEXITXX2-l.DO)/G5) 

AMU1=DARSINC1 . DO/MEXIT) 

PSI1=PAI2+AMU1 

ZETAF=G5XDATAN(DSQRT(MFIN**2-1 .D0)/G5) 

PSIF=PSI1+ZETA1-ZETAF 

NU0=0. 

TETLIM=NUFUNC(MEXIT)+PAI2-NUO 

PSILIM=TETLIM 

TETSYM=PAI-2 . DOXTETLIM 

GOREM=l . D0+G1XMEXITXX2 

RH01=RHOO/GOREMXXG6 

Vl=MEXITXAO/DSQRT(GOREM) 

Pl=PO/GOREM*x(G/(G-l .DO)) 

Tl=TO/DSQRT(GOREM) 

YYC=2.D0XPAIXYC 
IFCDELTA.EQ.O.) YYC=l.DO 
MD0T1=YYCXXCXRH01XV1 

PRINT 21,AR,AV,AN,G,RH00,N0,P0,T0,A0,D 
F0RMAT(/1X, ’THERMODYNAMIC DATA: '/ 

1 IX, ’AR,AV,AW,G=’,2X,2D14.5,2F9.3/ 

2 IX, ’RHOO,NO,PO,TO,AO,D=’,6D13.5) 

PRINT 22,XC,YC,MEXIT,RH01,P1,T1,V1,MD0T1,PSI1XDEG,PSIFXDEG, 

1 PSILIMXDEG 

F0RMAT(/1X, 'CORNER DATA: XC,YC= ’ , 2F9 .2/ 

1 IX, 'EXIT CONDITIONS: ’, 

2 2X, ’MEXIT,RH01,P1,T1,V1,MD0T1=’,F9.3,5D13.4/ 

3 IX, ’CENTERED FAN LIMITS:’, 

4 2X, ’PSI1,PSIF,PSILIM=’,3F10.3) 

PRINT 23, DELTA, KFO,JMAX,ITERO,DYO,YMAX, STAB, DDY 

F0RMAT(/1X, ’INTEGRATION DATA. SYMMETRY INDEX: DELTA= ’ , F4 . 1/ 



1 

2 

3 

5 

6 



RETURN 

END 



IX, ’NUMBER OF POINTS IN X AND Y DIRECTIONS: KFO,JMAX=' 

215/ 

IX, ’MAX. NUM. OF ITERATIONS ITER0=’,I5/ 

IX, ’INITIAL Y-STEP AND MAXIMUM Y: DYO ,YMAX= ’ , 2D14 .5/ 

1X,’Y-STEP STABILITY FACTORS STAB, DDY= ’ , 2F7 . 3) 



ABC IN 



SUBKUUliNE BtUiN 

SUBROUTINE NUMBER 2 

IMPLICIT REALX8(A-H,L-Z,$) 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MNa01),MUN(101), 

3 TETAN(101),BN(101),XTEMP(101) 
COMMON/THICKY/XTH(1002),TH(1002) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI,PAI2, DEG,XC, YC, MEXIT, MFIN, YMAX, DYO , DY, DYNEXT , 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 
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C 
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c 



COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

1 KF2,IDEL, JDEL,JYXI, JXI,ILEAD,ILEADF,KCLEAD 

COMMON /RON/YF,YN,DXF,DXN 

COMMON /CHARAC/XCHARFC 92 ) , YCHARF( 92) , XCHARNC 92) , YCHARN( 92 ) , 

1 RMCARFC 92) , RPCARF( 92) , RMCARN( 92) , RPCARN( 92) , 

2 TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 

4 MCHARK92) 

COMMON /ICHARA/KCHARP, KCHARM, KCHARO 

DEFINE INITIAL CHARACTERISTIC PARAMETERS. USE INTERPOLATION OF 
RIEMANN INVARIANT ACROSS THE FAN, 

KCHARP=IDINT(FRACGxDFL0AT(KF0-l)+l.D-6)+l 
KCHAR0=KCHARP+1 
KCHARM=KCHARO-KCHARP 
IF(KCHARP.LT.2 ) CALL FINC200) 

IF(KCHAR0.GT.92) CALL FIN(210) 

IF(KCHARM,LT. 1) CALL FIN(205) 

NU1=NUFUNC(MEXIT) 

RMl=NUO 
TET1=RM1-NU1 
RP0=NU1-TET1 
NUFIN=NUFUNC(MFIN) 

RPFIN=NUFIN-(RM1-NUFIN) 

DRP=(RPFIN-RPO)/DFLOAT(KCHARP-l) 

DO 1 KC=1,KCHARP 
RPl=RPO+DRP*DFLOAT(KC-l) 

CALL RFUNC(RM1,RP1,M1,MU1,TETA1) 

MCHARI(KC)=M1 
CONTINUE 

DATA FOR C+ CHARACTERISTICS. 

THE RIEMANN INVARIANTS ARE DEFINED IN SUCH A WAY THAT BOTH VANISH 
INFINITE MACH NUMBER. 

RMl=NUO 

DO 2 KC=1,KCHARP 
CSIGNF(KC)=l.DO 
XCHARF(KC)=XC 
YCHARF(KC)=YC 

IF(MCHARKKC) .EQ.O. ) CALL FIN(231) 

NU=NUFUNC(MCHARI(KC) ) 

TET=RM1-NU 

RP1=NU-TET 

CALL RFUNC(RM1,RP1,M1,MU1,TETA1) 

MCHARF(KC)=M1 
MUCARF(KC)=MU1 
TCHARF(KC)=TETA1 
RMCARF(KC)=RM1 
RPCARF(KC)=RP1 
2 CONTINUE 

: DATA FOR C- CHARACTERISTICS. 

KC1=KCHARP+1 
XCHARF(KC1)=0.8D0XXC 
DO 21 KC=KC1, KCHARO 
CSIGNF(KC)=-1 .DO 
MCHARI(KC)=MEXIT 

MUCARF(KC)=DARSIN(1.DO/MCHARI(KO) 

TCHARF(KC)=PAI2 

YCHARF(KC)=YC 

MCHARF(KC)=MEXIT 

RMCARF(KC)=RM1 

RPCARF(KC)=NUFUNC(MEXIT)-(TCHARF(KC)-TETLIM) 

21 CONTINUE 

: DEFINE GRID AND INITIAL CONDITIONS AT EXIT PLANE. 

KFAN=KCHARP-1 

ILEAD=KFO-KFAN 

KCLEAD=1 

KF=KFO 

KF2=2XKF 

YF=YC 

DO 3 1=1, KF 
KC=KCLEAD+I-ILEAD 
IF(KC.GT.KCHARP) CALL FIN(241) 
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PRINT' 



IF(KC.GE.l) GO TO 31 
XF(I) = DFLOAT(I-l)XXC/DFLOAT< ILEA D-D 
MF(I)=MEXIT 
TETAF(I)=PAI2 
GO TO 32 

31 CONTINUE 
XFf I ) =XC 
MF(D=MCHARF(KC) 

TETAF(I)=TCHARF(KC) 

32 CONTINUE 

RMF(I)=NUFUNC(MF(D)+(TETAF(I)-TETLIM) 
RPF(I)=NUFUNC(MF(D)-(TETAF<I)-TETLIM) 

MUF(D = DARSIN(1.D0/MFCD) 

BF(I)=0. 

3 CONTINUE 
DY=DY0 

DO A KC=1,KCHAR0 
CSIGNN(KC)=CSIGNF(KC) 

4 CONTINUE 
DO 5 1=1, KN 
BN(I)=0. 

5 CONTINUE 
RETURN 
END 

i5UBKUUIlNh"PKlNI 

: SUBROUTINE NUMBER 3 

IMPLICIT REALX8(A-H,L-Z,$) 

REAL5(4 XPL,YPL 

COMMON /PLUME/XPL ( 1 002, 10),YPL( 1002) 

COMMON /IPLUME/KPL,ITYPL(10) 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3 TETAN(101),BN(101),XTEMPa01) 

COMMON/THICKY/XTH(l 002), TH( 1002) 

REALX4 YXI,XI,XIPM,XIGRP,XIAPP,XIF 

COMMON /THICKX/YXI(20),XI(101,20),XIPM(101,20),XIGRP(101,20) 

1 ,XIAPP(101,20),XIF(101,20) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,611,G12,G13,G14, 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI , PAI2, DEG, XC, YC, MEXIT, MFIN, YMAX, DYO, DY, DYNEXT, 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 

COMMON /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), 

1 RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), 

2 TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 

4 MCHARK92) 

COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

1 KF2,IDEL, JDEL,JYXI, JXI,ILEAD,ILEADF,KCLEAD 

COMMON /ROW/YF,YN,DXF,DXN 

SUM=0. 

KF1=KF-1 

DO 10 1=1, KFl 

DX=XF(I+1)-XF(I) 

GOREM=l .D0+G1XMF(I)XX2 
GOREM1=1.DO+G1XMF(I+1)XX2 

RATEM = RHOOXAO?(MF(I )XDSIN(TETAF( I ) )/GOREM XX(G6 + 0.5D0) 

RATEP = RHO05(A03fMF(I + l)XDSIN(TETAF(I + l))/GOREMlxx(G6 + 0.5D0) 
SUM=SUM+DXX(RATEM+RATEP)/2 . DO 

10 CONTINUE 
YYF=2.D0XPAIXYF 
IF(DELTA.EQ.O. ) YYF=1.D0 
MDOTFR=YYFXSUM/MDOTl 

PRINT 11, J,KCLEAD,KF,ILEAD,YF,DY,XF(KF),MF(KF),MDOTFR 

11 FORMATCIX, ’J,KCLEAD,KF,ILEAD,YF,DY,XF(KF),MBOUND,MDOTR=', 

1 4I5,5D12.4) 



C 

C 

C 



PRINT FLOH FIELD AT Y=YF 
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IF(J .EQ. JMAX) JXI=MINO(JXI, JYXI) 

IFCJ.EQ.l .OR. J.EQ.JMAX) GO TO 121 
IF(JXI.GT.JYXI) GO TO 120 
IF(YXI(JXI).GT.YF+0.5D0*DY) GO TO 120 
121 CONTINUE 

YXI(JXI)=YF 
CALL OPACX 

C COMPUTE MACH NUMBER FOR CYLINDRICAL EXPANSION MCYL . 
F=(YF/YC)x(G7x(l .D0+G1XMEXIT**2))*3(G2/MEXIT 
CALL AREAF(F,MCYL) 

PRINT 22,JXI,KCLEAD,ILEAD,KF,MCYL,YF 
22 FORMATC/IX, 'PRINTING NUMBER JXI , KCLEAD, ILEAD, KF= ' , 414, 
1 5X, *MCYL,YF=',2D14.5/) 

PRINT 1 



FORMATC/IX, ' I ',' KC 


',' XFCI) ',' 


TETAFCI) 




1 


' MFCI) ',' 


MAB 


ff 


2 


' MABI ',' 


MOBI 


ff 


3 


' XICI) ',' 


XIGRPCI) 


ff 


4 


' XIAPPCI) ',' 


XIPMCI) 


'/) 


IDEL1=IDEL 

IFCJ.EQ.l. OR. J.EQ.JMAX) 


IDEL1=1 







DO 20 I=1,KF,IDEL1 
KC=KCLEAD+(I-ILEAD) 

IF(KC.LT. KCLEAD) KC=0 

M0B=1.D10 

M0BI=1 .DIO 

MAB = 1 .DIO 

MABI=1 .DIO 

MPM=MF(I) 

IF(KC.EQ.O) GO TO 23 
MOB=MCHARI(KC) 

IFCJ.EQ.l) GO TO 23 

PSIPM=PAI2-DATAN((XF(I)-XC)/(YF-YO) 

ZETA=PSI1+ZETA1-PSIPM 

MPM=DSQRT( (G55<DTAN{ZETA/G5) )XX2+1 . DO) 

CALL MATCH(I,MOB,MAB,MOBI,MABI) 

23 CONTINUE 

PRINT 21,I,KC,XF(I),TETAF(I)3(DEG,MF(I),MAB,MABI,M0BI, 

1 XKI, JXI),XIGRP(I,JXI),XIAPP(I,JXI),XIPM(I,JXI) 

21 FORMAT(1X,2I4,10D12.4) 

20 CONTINUE 

IFCJ.EQ.l) GO TO 120 
IFCJ.EQ.JMAX) GO TO 120 
JXI=JXI+1 
120 CONTINUE 

IFC J .LT. JMAX) GO TO 200 
PRINT 101 

101 FORMATC'l') 

PRINT 102 

102 FORMATC IX, 'RADIAL MOLECULAR THICKNESS J , XTHC J ) , THC J ) = '/ ) 
PRINT 202, CJ J, XTHC JJ), THC JJ),JJ=1, JMAX) 

202 FORMATC/5CI5,D11.4,D10.3)) 

PRINT 101 

PRINT 103, CIPL,ITYPLCIPL),IPL=1,KPL) 

103 FORMATCIX, 'PLUME TYPES IPL , ITYPL C IPL ) = ' , 

1 2C/1X,5C5X,2I4))) 

PRINT 104 

104 FORMATCIX, 'PLUME POINTS J ,YPL CJ ) ,XPL C J , 1 ), XPL C J , 2) ,...='/ ) 
JDEL1=1 

DO 203 JJ=1, JMAX, JDELl 

PRINT 204, JJ,YPLCJJ),CXPLCJJ, IPL),IPL=1,KPL) 

204 FORMATC1X,I5,2X,E12.4,10E11 .3) 

203 CONTINUE 

C WRITE ON TAPE? FOR SUBSEQUENT PLOTTING. 

C NO MORE THAN 80 CHARACTERS PER LINE ON TAPE?. 

WRITEC?,205) JMAX,KPL 

205 FORMATC8I10/8I10) 

WRITEC?,205) C ITYPL C IPL), IPL=1,KPL) 

DO 210 JJ=1,JMAX 

WRITEC?,211) YPLCJJ),CXPLCJJ,IPL),IPL=1,KPL) 

211 F0RMATC6E13.6/2X,6E13.6/2X,6E13.6/2X,6E13.6) 

210 CONTINUE 
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JXI0,KF0=',2I8) 



: WRITE LATERAL (X) OPACITIES 

JXI0=JXI 

WRITE(7,205) JXI0,KF0 
PRINT 226, JXI0,KF0 

226 F0RMAT(///1X, 'LATERAL (X) OPACITIES 
DO 220 JXI=1,JXI0 

WRITE(7,221) JXI,YXI(JXI) 

221 FORMAT(I10,E13.6) 

PRINT 227, JXI,YXI(JXI) 

227 F0RMAT(//1X, • JXI,YXI( JXI)=', 18, E15.6/) 

DO 225 1=1, KFO 

WRITE(7,211) XIF(I,JXI),XI(I,JXI),XIPM(I,JXI),XIGRP(I,JXI), 
XIAPPd, JXI) 



PRINT 211, XIFCI, JXI),Xia,JXI),XIPM(I, JXI),XIGRPa,JXI), 



225 

220 

200 



CONTINUE 

CONTINUE 

CONTINUE 

RETURN 



XIAPPd, JXI) 



f=lN 
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SUBROUTINF'PINnFIN) 
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c 


SUBROUTINE NUMBER 4 




JET0454 


c 


STOP WHEN ERROR IS DETECTED. 




JET0455 




IMPLICIT REALX8(A-H,L-Z,$) 




JET0456 




PRINT 1,IFIN 




JET0457 


1 


F0RMAT(/1X, 'FIN CODE IFIN=',I6/) 




JET0458 


c 


INDUCE ERROR IN ORDER TO GENERATE TRACING 


OF CALIINO SUBROUTINES. 


JET0459 




X=-1.D0 




JET0460 




Y=X+DSQRT(X) 




JET0461 




IF(IFIN.LE.O) GO TO 100 




JET0462 




STOP 
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1 


00 RETURN 
END 
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SUBROUTINE NUMBER 5 JET0467 

IMPLICIT REAL5«8(A-H,L-Z,$) JET0468 

COMMON /VECS/XFd01),RMFd01),RPFd01),MFd01),MUFd01), JET0469 

1 TETAF(101),BFd01), JET0470 

2 XNd01),RMN(101),RPN(101),MNd01),MUNd01), JET0471 

3 TETANd01),BNd01),XTEMPd01) JET0472 

C0MM0N/THICKY/XTHd002),THd002) JET047 3 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, JET0474 

1 G16,G17,G18,G19,G20 JET0475 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, JET0476 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NU0, JET0477 

2 TETSYM,TETLIM,DDY,DYMAX JET0478 

COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 JET0479 

COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, JET0480 

1 KF2,IDEL, JDEL, JYXI, JXI,ILEAD,ILEADF,KCLEAD JET0481 

COMMON /ROW/YF,YN,DXF,DXN JET0482 

COMMON /CHARAC/XCHARF( 92 ) , YCHARF( 92 ) , XCHARNC 92) , YCHARNC 92) , JET0483 

1 RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), JET0484 

2 TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), JET0485 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), JET0486 

4 MCHARK92) JET0487 

COMMON /ICHARA/KCHARP,KCHARM,KCHARO JET0488 

JET0489 

ADVANCE FLOW FIELD FROM YF TO YN JET0490 

IM=KF JET0491 

IP=KF JET0492 

YN=YF+DY JET0493 

KN=KF0 JET0494 

SEMI-INVERSE INTEGRATION FOR FAN POINTS. JET0495 

CALL SEMINV JET0496 

NEW GRID POINTS (JUST INVERSE MARCHING). JET0497 

CALL GRIDN JET0498 

LOAD FLOW VARIABLES FROM SEMI-INVERSE INTEGRATION INTO VECTORS JET0499 

CALL LOADC JET0500 

CHARACTERISTIC SCHEME INTEGRATION FOR INNER POINTS (INVERSE MARCH). JET0501 
CALL INVMAR JET0502 

RETURN JET0503 

END JET0504 



C 

c 
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C 

c 
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SUBROUTINE INVMAR 
SUBROUTINE NUMBER 6 

IMPLICIT REALX8(A-H,L-Z,$) 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3 TETAN(101),BN(101),XTEMP(101) 

COMMON/THICKY/XTH( 1002), TH( 1002) 

COMMON /GAMA/6, G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1 616,G17,G18,G19,G20 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 
COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

1 KF2,IDEL, JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON /ROW/YF,YN,DXF,DXN 

INTEGRATION WITH INVERSE CHARACTERISTICS FOR NEW P0INT(X4,Y4) . 

OLD POINTS ARE (X1,Y1), (X2,Y2) . 

XI IS OBTAINED BY INVERSE C- FROM X4 
X2 IS OBTAINED BY INVERSE C+ FROM X4 
NOTE THAT XI MAY BE NEGATIVE (E. G. WHEN X4=0). 

KN1=ILEAD-1 

IF(KNl.LE.O) CALL FIN(601) 

DO 1000 1=1, KNl 
14 = 1 

X4=XN(I) 

Y4=YN 

IF4=(IM+IP)/2 

CALL INTERP(0,IF4,KF,X4,XF,RM4,RMF,RP4,RPF) 

CALL RFUNC(RM4,RP4,M4,MU4,TETA4) 

M14=M4 

MU14=MU4 

TETA14=TETA4 

M24=M4 

MU24=MU4 

TETA24=TETA4 

Y1=YF 

Y2=YF 

Y14=(Y1+Y4)/2.D0 
Y24=(Y2+Y4)/2.D0 
X1=1.D10 
X2=1.D10 
RM4=1 .DIO 
RP4=1 .DIO 
ITER=0 
GO TO 2 

CORRECTOR 

ITER=ITER+1 

AVERAGED PROPERTIES ON C-( 14 ) , C+( 24) CHARACTERISTICS. 
RM14=(RM1+RM4)/2.D0 
RP14=(RP1+RP4)/2.D0 
RM24=(RM2+RM4)/2.D0 
RP24=(RP2+RP4)/2.D0 

M14,MU14,TETA14, M24 , MU24 . TETA24 AVERAGED ON C-,C+ CHARACTERISTICS. 
CALL RFUNC(RM14,RP14,M14,MU14,TETA14) 

CALL RFUNC(RM24,RP24,M24,MU24,TETA24) 

CONTINUE 
NEW XI, X2 
X10=X1 
X20=X2 

X1=X4-DY/DTAN(TETA14-MU14) 

X2=X4-DY/DTAN(TETA24+MU24) 

IF(X2.LT.O.) CALL FIN(670) 

D14=DSQRT((X1-X4)XX2+DYXX2) 

D24=DSQRT( (X2-X4)xx2+DY*x2) 

INTERPOLATE OLD DISTRIBUTION FOR RM1,RP1, RM2,RP2 AT XI, X2. 

CALL INTERP(0,IM,KF,X1,XF,RM1,RMF,RP1,RPF) 

CALL INTERP(0,IP,KF,X2,XF,RM2,RMF,RP2,RPF) 



JET0505 
JET0506 
JET0507 
JET0508 
JET0509 
JET0510 
JET0511 
JET0512 
, JET0513 
JET0514 
JET0515 
JET0516 
JET0517 
JET0518 
JET0519 
JET0520 
JET0521 
JET0522 
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JET0526 
JET0527 
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JET0531 
JET0532 
JET0533 
JET0534 
JET0535 
JET0536 
JET0537 
JET0538 
JET0539 
JET0540 
JET0541 
JET0542 
JET0543 
JET0544 
JET0545 
JET0546 
JET0547 
JET0548 
JET0549 
JET0550 
JET0551 
JET0552 
JET0553 
JET0554 
JET0555 
JET0556 
JET0557 
JET0558 
JET0559 
JET0560 
JET0561 
JET0562 
JET0565 
JET0564 
JET0565 
JET0566 
JET0567 
JET0568 
JET0569 
JET0570 
JET0571 
JET0572 
JET0573 
JET0574 
JET0575 
JET0576 
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NO NEED FOR RE-AVERAGING SINCE IT INTRODUCES ONLY HIGHER ORDER 
CHANGES INTO THE ITERATION SCHEME. 

INTEGRATE THE CHARACTERISTIC EQUATIONS FOR RMA,RPA AT XA,YA. 

RMAO=RMA 

RPAO=RPA 

RMA=RM1+DELTA3<DSIN(TETA1A)*D1A/(M1AXY1A) 

RPA=RP2+DELTA*DSIN(TETA2A)XD2A/(M24*Y24) 

CONVERGENCE TEST 

EPS=(DABS(X1-X10)+DABS(X2-X20))/DY+DABS(RMA-RM40)+DABS(RP4-RPA0) 

IF(ITER.GT.ITERO) GO TO 10 

IF(EPS.GT.EPSIL) GO TO 1 

RMN(I)=RM4 

RPN(I)=RP4 

CALL RFUNC(RM4,RP4,MN(I),MUN(I),TETAN(D) 

1000 CONTINUE 
RETURN 
CONTINUE 

PRINT 11,I4,KN,IF4,IM,IP,KF,ITER,ITER0,EPS,EPSIL,X1,X2,X4,M14,M24 
FORMATdX, 'SUBR. INVMAR. 14, KN, IF4, IM, IP, KF, ITER, ITER0= • ,815/ 

IX, 'EPS,EPSIL,X1,X2,X4,M14,M24=*,7D14.6/) 



10 

11 



CALL FIN(611) 
RETURN 
END 

SUBROUTINE 5EMINV 



SSMNi/' 
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SUBROUTINE NUMBER 7 

IMPLICIT REAL*8(A-H,L-Z,$) 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3 TETAN(101),BN(101),XTEMP(101) 
COMMON/THICKY/XTH(1002),TH(1002) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 
COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

1 KF2,IDEL, JDEL, JYXI, JXI,ILEAD,ILEADF,KCLEAD 

COMMON /ROW/YF,YN,DXF,DXN 

COMMON /CHARAC/XCHARF( 92) , YCHARFC 92) , XCHARNC 92 ) , YCHARN( 92) , 

1 RMCARF(92) , RPCARFC 92) , RMCARN( 92) , RPCARN( 92) , 

2 TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 

4 MCHARK92) 

COMMON /ICHARA/KCHARP, KCHARM, KCHARO 

COMPUTE NEW POINT (X4,Y4), BY PASSING A C+ CHARACTERISTIC 
THROUGH OLD POINT (X2,Y2). BOTH POINTS ARE ON CHARACTERISTIC LINE 
NUMBER KC. 

IM = 1 

DO 100 KC=1, KCHARO 
IF(CSIGNN(KC) .EQ.O.) GO TO 100 

PREDICTOR 

Y1=YF 

Y2=YF 

Y4=YN 

Y14=(Y1+Y4)/2.D0 

Y24=(Y2+Y4)/2.D0 

X2=XCHARF(KC) 

RM2=RMCARF(KC) 

RP2=RPCARF(KC) 

M2=MCHARF(KC) 

MU2=MUCARF(KC) 

TETA2=TCHARF(KC) 

M14=M2 

MU14=MU2 

TETA14=TETA2 

M24=M2 

MU24=MU2 
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TETA24=TETA2 
X4=1.D10 
X1=1.D10 
RM4=1.D10 
RP4=1.D10 
ITER=0 
GO TO 2 
C 

C CORRECTOR 
C 

1 ITER=ITER+1 

C AVERAGED PROPERTIES ON C-(14),C+(24) CHARACTERISTICS. 
RM14=(RM1+RM4)/2.D0 
RP14=(RP1+RP4)/2.D0 
RM24=(RM2+RM4)/2.D0 
RP24=(RP2+RP4)/2.D0 

C M14,MU14,TETA14, M24,MU24,TETA24 AVERAGED ON C-,C+ CHARACTERISTICS. 
CALL RFUNC(RM14,RP14,M14,MU14,TETA14) 

CALL RFUNC(RM24,RP24,M24,MU24,TETA24) 

2 CONTINUE 
C NEW X4,X1 

X40=X4 

X10=X1 

X4=X2+DY/DTAN(TETA24+CSIGNF(KC)*MU24) 

X1=X4-DY/DTAN(TETA14-CSIGNF(KC)3(MU14) 

D14=DSQRT((X1-X4)XX2+DYXX2) 

D24=DSQRT( (X2-X4)x*2+DYxx2) 

C INTERPOLATE OLD DISTRIBUTION FOR RM1,RP1, AT XI. 

CALL INTERP(0,IM,KF,X1,XF,RM1,RMF,RP1,RPF) 

IF(J.GT.l) GO TO 22 
IF(CSIGNF(KC) .LT.O.) GO TO 22 
RP1=RP2 
22 CONTINUE 

C NO NEED FOR RE-AVERAGING SINCE IT INTRODUCES ONLY HIGHER ORDER 
C CHANGES INTO THE ITERATION SCHEME. 

C INTEGRATE THE CHARACTERISTIC EQUATIONS FOR RM4,RP4 AT X4,Y4. 
RM40=RM4 
RP40=RP4 

IF(CSIGNF(KC) .LT.O. ) GO TO 21 

RM4 = RM1+DELTAXDSIN(TETA14)5(D14/(M14*Y14) 

RP4 = RP2+DELTA*DSIN(TETA24)*D24/(M243(Y24) 

GO TO 20 
21 CONTINUE 

RM4=RM2+DELTAXDSIN(TETA24)XD24/(M24*Y24) 

RP4 = RP1 + DELTA5«DSIN(TETA14)5«D14/(M14*Y14) 

20 CONTINUE 
C CONVERGENCE TEST 

EPS=(DABS(X4-X40)+DABS(X1-X10))/DY+DABS(RM4-RM40)+DABS(RP4-RP40) 
IF(ITER.GT.ITERO) GO TO 10 
IF(EPS.GT.EPSIL) GO TO 1 
CSIGNN(KC)=CSIGNF(KC) 

IFCX4.GT.0.) GO TO 30 
RMSAVE=RM4 
RM4=RP4+TETSYM 
RP4=RM4-TETSYM 
CSIGNN(KC)=-1 .DO 
30 CONTINUE 

RMCARN(KC)=RM4 

RPCARN(KC)=RP4 

CALL RFUNC(RM4,RP4,M4,MU4,TETA4) 

TCHARN(KC)=TETA4 

XCHARN(KC)=DABS(X4) 

YCHARN(KC)=Y4 
MUCARN(KC)=MU4 
MCHARN(KC)=M4 
100 CONTINUE 
RETURN 

10 CONTINUE 

PRINT 11,KC,KCHAR0,IM,KF, ITER, ITERO , EPS, EPSIL , XI , X2, X4 , M14 ,M24 

11 FORMATdX, 'SUBR. SEMINV. KC, KCHARO , IM, KF, ITER, ITER0= ', 615/ 

1 IX, ’EPS,EPSIL,X1,X2,X4,M14,M24=',7D14.6/) 

CALL FIN(711) 



JET0649 

JET0650 

JET0651 

JET0652 

JET0653 

JET0654 

JET0655 

JET0656 

JET0657 

JET0658 

JET0659 

JET0660 

JET0661 

JET0662 

JET0663 

JET0664 

JET0665 

JET0666 

JET0667 

JET0668 

JET0669 

JET0670 

JET0671 

JET0672 

JET0673 

JET0674 

JET0675 

JET0676 

JET0677 

JET0678 

JET0679 

JET0680 

JET0681 

JET0682 

JET0683 

JET0684 

JET0685 

JET0686 

JET0687 

JET0688 

JET0689 

JET0690 

JET0691 

JET0692 

JET0693 

JET0694 

JET0695 

JET0696 

JET0697 

JET0698 

JET0699 

JET0700 

JET0701 

JET0702 

JET0703 

JET0704 

JET0705 

JET0706 

JET0707 

JET0708 

JET0709 

JET0710 

JET0711 

JET0712 

JET0713 

JET0714 

JET0715 

JET0716 

JET0717 

JET0718 

JET0719 

JET0720 



29 



JETPR 



FORTRAN A1 



RETURN 
END 



RFUHC 



JET0721 

■IFT.nT2?.. 



SUBROUTINE RFUNC(RM,RP,M,MU,TETA) JET0723 

C SUBROUTINE NUMBER 8 JET072A 

IMPLICIT REALX8(A-H,L-Z,$) JET0725 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), JET0726 

1 TETAF(101),BF(101), JET0727 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), JET0728 

3 TETAN(101),BN(101),XTEMP(101) JET0729 

COMMON/THICKY/XTH( 1002 ),TH( 1002) JET0730 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,611,G12,G13,G1A,G15,JET0731 

1 G16,G17,G18,G19,G20 JET0732 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, JET0733 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NU0, JET0734 

2 TETSYM,TETLIM,DDY,DYMAX JET0735 

COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 JET0736 

COMMON /IPAR/JMAX,KF0,ITER0,KF,KN,IM,IP,J, JET0737 

1 KF2,IDEL, JDEL, JYXI,JXI,ILEAD,ILEADF,KCLEAD JET0738 

COMMON /ROW/YF,YN,DXF,DXN JET0739 

C JET0740 

C COMPUTE M,MU,TETA AT A POINT, AS FUNCTION OF RIEMANN INVAR. RM,RP. JET0741 

TETA=(RM-RP)/2.D0+TETLIM JET0742 

NU =(RM+RP)/2.D0 JET0743 

C NU=NU0-(G5XARCTAN(G5XQ)-ARCTAN(Q)), where Q=(MXX2-1)xx(-1/2) JET0744 

C FIND Q(NU), AND HENCE M(NU), THROUGH NEWTON RAPHSON ITERATIONS. JET0745 

Q=-(NU-NU0)/(G4-1.D0) JET0746 

IFCQ.LE.O.) CALL FIN(801) JET0747 

ITER=0 JET0748 

1 ITER=ITER+1 JET0749 

QF=Q JET0750 

DNUDT=-(G4-1 . D0)/( ( 1 . D0+G4XQXX2)X( 1 . D0+QXX2) ) JET0751 

DNU=NU-(NU0-(G5XDATAN(G5*Q)-DATAN(Q))) JET0752 

Q=Q+DNU/DNUDT JET0753 

IFCQ.LE.O.) CALL FINCSll) JET0754 

EPS=DABS(Q-QF)/Q JET0755 

IFCITER.GT.ITERO) GO TO 10 JET0756 

IFCEPS.GT.EPSILXl.D-3) GO TO 1 JET0757 

M=DSQRT(1 .DO+1 .D0/QXX2) JET0758 

MU=DARSIN(1 .DO/M) JET0759 

RETURN JET0760 

10 CONTINUE JET0761 

CALL FINC810) JET0762 

RETURN UJTFQ P JET0763 

END ' JET0764 
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SUBROUTINE INTERPC JNEW, I,KGRID,X,XVEC,RM,RMVEC,RP,RPVEC) 

SUBROUTINE NUMBER 9 

IMPLICIT REALX8(A-H,L-Z,#) 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3 TETAN(101),BN(101),XTEMP(101) 

COMMON/THICKY/XTHCl 002 ),TH( 1002) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1 G16 , G17 , G18 , G1 9 , G20 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1 STAB, DELTA, PSIl , PSIF, ZETAl , SIGMA, FRACG, EPSIL , NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RHOO , NO , PO , TO , AO , MDOTl 
COMMON /IPAR/JMAX,KFO,ITERO,KF,KN, IM,IP, J, 

1 KF2,IDEL, JDEL,JYXI, JXI,ILEAD,ILEADF,KCLEAD 

COMMON /ROW/YF,YN,DXF,DXN 
DIMENSION XVECC 1 ) , RMVECC 1 ) , RPVECC 1 ) 



FIND I SUCH THAT XVECC I ). LE . X . AND . XVECC I+l ) .GE .X 
FIND RM,RP BY LINEAR INTERPOLATION. 

NOTE THAT X MAY BE NEGATIVE. 

IFCDABSCX).LE.XVECCKGRID)) GO TO 901 
PRINT 900,X,KGRID,XVECCKGRID) 
F0RMATC/1X,D15.7, I10,4X,D15.7/) 

CALL FINC900) 

901 CONTINUE 

KG2=2XKGRID 
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I0=MIN0(I,KGRID-2) 

ICOUNT=0 
1 1 = 10 

SIGN1=1.D0 
IF(I.GE.l) GO TO 10 
SIGN1=-1.D0 
I=-I+2 

10 CONTINUE 

IF(I.GT.KGRID) CALL FIN(901) 
XXl=SIGNlxxVECm 
11 = 1 

IF(XXl.LE.X) GO TO 11 
10 = 10-1 

IC0UNT=IC0UNT+1 

IF(IC0UNT.GT.KG2) CALL FIN(911) 

GO TO 1 

11 CONTINUE 
1 = 10+1 
SIGN2=1.D0 
IF(I.GE.l) GO TO 12 
SIGN2=-1.D0 
I=-I+2 

12 CONTINUE 

IF(I.GT.KGRID) CALL FIN(912) 
XX2=SIGN2XXVEC(I) 

12=1 

IF(XX2.GE.X) GO TO 13 
10 = 10+1 

IC0UNT=IC0UNT+1 

IF(IC0UNT.GT.KG2) CALL FIN(913) 

GO TO 1 

13 CONTINUE 
F1=(XX2-X)/(XX2-XX1) 

F2=1.D0-F1 

IF(F1.LT.0.) CALL FIN(991) 
IF(F2.LT.O.) CALL FIN(992) 
RM1=RMF(I1) 

RP1=RPF(I1) 

RM2=RMF(I2) 

RP2=RPF(I2) 

IF( SIGNl . LT. 0 . ) RM1=RPF( I1)+TETSYM 
IFCSIGNl.LT.O. ) RP1=RMF(I1)-TETSYM 
IFCSIGN2.LT.0. ) RM2=RPF( I2)+TETSYM 
IFCSIGN2.LT.0.) RP2=RMF( I2)-TETSYM 
RM=F15(RM1 + F25(RM2 
RP = F15(RP1+F2XRP2 
RETURN 
END 






JET0793 

JET079A 

JET0795 

JET0796 

JET0797 

JET0798 

JET0799 

JET0800 

JET0801 

JET0802 

JET0803 

JET080A 

JET0805 

JET0806 

JET0807 

JET0808 

JET0809 

JET0810 

JET0811 

JET0812 

JET0813 

JET081A 

JET0815 

JET0816 

JET0817 

JET0818 

JET0819 

JET0820 

JET0821 

JET0822 

JET0823 

JET0824 

JET0825 

JET0826 

JET0827 

JET0828 

JET0829 

JET0830 

JET0831 

JET0832 

JET0833 

JET083A 

JET0835 

JET0836 

JET0837 

JET0838 

JET0839 

JET0840 



C 

C 



3mTRUUIlNE"TNrtRXLJNtW,ii',VAR0';’V7nT7TCUKiU;X07XVh’C'J JFHJW 

SUBROUTINE NUMBER 10 JET08A2 

IMPLICIT REALX8(A-H,L-Z,$) JET0843 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), JET0844 

1 TETAF(101),BF(101), JET0845 

2 XN(101),RMN(101),RPN(101),MN(101),MUNa01), JET0846 

3 TETAN(101),BN(101),XTEMP(101) JET0847 

COMMON/THICKY/XTHC 1 002), TH( 1002) JET0848 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, JET0849 

1 G16,G17,G18,G19,G20 JET0850 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, JET0851 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO'> JET0852 

2 TETSYM,TETLIM,DDY,DYMAX JET0853 

COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 JET0854 

COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, JET0855 

1 KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD JET0856 

COMMON /ROW/YF,YN,DXF,DXN JET0857 

DIMENSION VAR(1),XVEC(1) JET0858 

FIND XO AND II SUCH THAT XVECCIl )<X0<XVEC( Il+l ) , AND XO CORRESPONDS JET0859 
TO THE LOCATION AT WHICH VARO IS A LINEAR INTERPOLATION OF VAR(I). JET0860 

X0=1.D23 JET0861 

IFIRST=I1 JET0862 

IF(Il.GT.O) GO TO 10 JET0863 

IFIRST=KGRID-IABS(Il)+2 JET0864 
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10 CONTINUE 

DO 1 II=IFIRST,KGRID 
I = II 

IF(Il.GT.O) GO TO 11 
I=KGRID-II+2 

11 CONTINUE 

IF(I.LE.O) CALL FIN(lOOl) 
IF(I.GT.KGRID) CALL FIN(1002) 
IF(I.EQ.l) GO TO 1 

IF((VAR(I)-VAR0)X(VAR(I-1)-VAR0) .GT 
IF(VARd) .EQ.VAR(I-D) GO TO 1 
F1 = (VAR(I)-VAR0)/(VAR(I)-VAR(I-D) 
F2=1.D0-F1 

IF(F1.LT.0.) CALL FIN(lOll) 
IF(F2.LT.O.) CALL FIN(1012) 

X0 = F15«XVEC(I-1)+F2XXVEC(I) 

11 = 1-1 
GO TO 2 

1 CONTINUE 

2 CONTINUE 
RETURN 
END 

SUBRUUIlNb BREAK 



0.) GO TO 1 
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SUBROUTINE NUMBER 11 

IMPLICIT REAL*8(A-H,L-Z,$) 

REAL MB, MX, MY 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF<101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN<101), 

3 TETAN(101),BN(101),XTEMP(101) 

COMMON/THICKY/XTHd 002), TH( 1002) 

COMMON /GAMA/G,G1,G2,G5,GA,G5,G6,G7,G8,G9,G10,G11,G12,G13,G1A 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 
COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM, IP, J, 

1 KF2,IDEL, JDEL, JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON /ROW/YF,YN,DXF,DXN 



COMPUTE THE BREAKDOWN PARAMETER AT d-1/2, K-1/2) . 
YB = 0.5D05«(YF+YN) 

DYY=DY 

IM=2 

DO 1 1=2, KN 
Xl = XNd-l) 

X2=XNd) 

nxx=X2-xi 

IF(X2.GT.XF(KF)) GO TO 2 

CALL INTERP(0,IM,KF,X1,XF,RM1,RMF,RP1,RPF) 

CALL INTERP(0,IM,KF,X2,XF,RM2,RMF,RP2,RPF) 

CALL RFUNC(RM1,RP1,M1,MU1,TETA1) 

CALL RFUNC(RM2,RP2,M2,MU2,TETA2) 
MX=0.5D0x((MNd)-MNd-l)) + (M2-Ml))/DXX 
MY=0.5DOX((MNd)-M2) + (MNd-l)-Ml))/DYY 
MB = 0.25DOX(MNd-l)+MNd)+Ml+M2) 

TETAB = 0.25D0X(TETANd-l)+TETANd)+TETAl+TETA2) 

G0REM=MBXX2X(1 . D0+G1XMBXX2)XX(G6-1 . DO) 

GRAD=MXXDC0S(TETAB)+MY5«DSIN(TETAB) 

B=G20XGOREMXGRAD 

GO TO 3 

B=1 .D22 

BNd) = B 

CONTINUE 

BNd)=BN(2) 

RETURN 
END 



STORE IN BN(I) 



OFAC'j 
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SUBKUUliNb UHACY 

SUBROUTINE NUMBER 12 

IMPLICIT REALX8(A-H,L-Z,$) 

COMMON /VECS/XFd01),RMFd01),RPFd01),MFd01),MUFd01), 
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C 

C 



C 

c 

c 

c 

c 

c 

c 

c 

c 

c 



1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3 TETAN(101),BN(101),XTEMP(101) 

COMMON/THICKY/XTH( 1002), TH( 1002) 

COMMON /GAMA/G,G1,G2,G3,GA,G5,G6,G7,G8,G9,G10,G11,G12,G13,G1A,G15, 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI , PAI2, DEG, XC, YC, MEXIT, MFIN, YMAX, DYO, DY, DYNEXT, 

1 STAB, DELTA, PSIl, PSIF, ZETAl, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 
COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

1 KF2,IDEL,JDEL, JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON /ROW/YF,YN,DXF,DXN 

COMPUTE THE MOLECULAR THICKNESS AT END POINTS OF EACH ROM. 

IM=2 

XTH(J)=XF(KF) 

TH(J)=0. 

DTHO=NOXSIGMAXDY 
IF(J.EQ.l) GO TO 11 
J1=J-1 

DO 1 JJ=1,J1 
XX1=XTH(JJ) 

CALL INTERP(0,IM,KF,XX1,XF,RM1,RMF,RP1,RPF) 

CALL RFUNC(RM1,RP1,M1,MU1,TETA1) 

GOREM=l . D0+G1XM1XX2 

DTH=DTHO/GOREMXxG6 

TH(JJ)=TH(JJ)+DTH 

I CONTINUE 

II CONTINUE ^ 

return PLUM£S 

bUtfRUUTTTTE“POJHtS 

: SUBROUTINE NUMBER 13 

IMPLICIT REALX8(A-H,L-Z,$) 

REALXA XPL,YPL 

COMMON /PLUME/XPL(1002,10) ,YPL(1002) 

COMMON /IPLUME/KPL,ITYPL(10) 

DIMENSION VPLC92) 

COMMON /VECS/XFa01),RMF(101),RPF(101),MF(101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3 TETAN(101),BN(101),XTEMP(101) 

COMMON/THICKY/XTHC 1002), TH( 1002) 

REALX4 YXI,XI,XIPM,XIGRP,XIAPP,XIF 

COMMON /THICKX/YXI(20),XI(101,20),XIPM(101,20),XIGRP(101,20) 

1 ,XIAPP(101,20),XIFa01,20) 

COMMON /GAMA/G,G1,G2,G3,G<4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G1A,G15, 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI, PAI2, DEG, XC,YC, MEXIT, MFIN, YMAX, DYO, DY, DYNEXT, 

1 STAB, DELTA, PSIl, PSIF, ZETAl, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 
COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

1 KF2,IDEL,JDEL, JYXI, JXI,ILEAD,ILEADF,KCLEAD 

COMMON /ROH/YF,YN,DXF,DXN 

COMMON /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), 

1 RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), 

2 TCHARF(92),TCHARN(92),MUCARF(92) ,MUCARN(92), 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 

A MCHARK92) 

COMMON /ICHARA/KCHARP , KCHARM, KCHARO 
COMPUTE SPECIAL POINTS AT Y=YN, AND STORE THEM AS 
(XPL(J,IPL),YPL( J)=YN) . 

J IS THE MARCHING INDEX OF YN . 

IPL=1,2, . . . ,KPL IS THE "PLUME" INDEX. PRESENTLY KPL.LE.5 
VPL(IPL) IS A VALUE DEFINING THE "PLUME" CURVE. 

ITYPL(IPL) IS THE TYPE OF CURVE. IT DEFINES CURVES AS FOLLOWS: 

DO NOTHING 

REAL PLUME. IT IS THE BREAKDOWN SURFACE, DEFINED 
BY A CONSTANT VALUE OF THE BREAKDOWN PARAMETER B. 

SET VPL(IPL)=B. 



ITYPL(IPL)=0 

ITYPL(IPL)=1 
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C 

C 

C 

c 

c 

c 

c 

c 

c 

c 



ITYPL(IPL)=2 

ITYPL(IPL)=3 

ITYPL(IPL)=4 



ITYPL(IPL)=5 



CONSTANT MACH-NUMBER LINE. VPL(IPL)=M. 

A SINGLE STREAMLINE. VPL(IPL) IS SET TO THE EXIT 
X-COORDINATE OF THAT STREAMLINE. 

A SINGLE C+ CHARACTERISTIC LINE STARTING AT THE CORNER, 
VPL(IPL) IS SET TO THE INDEX KC OF THAT CHARACTERISTIC 
LINE. 

A CONSTANT LATERAL (X) OPACITY LINE. VPL(IPL) IS SET 
TO THE VALUE OF THE (CONSTANT) OPACITY. 



DEFINE ITYPL(IPL) AND VPL(IPL) 

KPL=10 

IF(KPL.GT.IO) CALL FINCISOI) 

DO 2000 IPL=1,KPL 

GO TO (2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), IPL 

2001 ITYPL(IPL)=4 
VPL(IPL)=1 
GO TO 2000 

2002 ITYPL(IPL)=4 
VPL(IPL)=KCHARP 
GO TO 2000 

2003 ITYPL(IPL)=4 
VPL(IPL)=19 
GO TO 2000 

2004 ITYPL(IPL)=4 
VPL(IPL)=31 
GO TO 2000 

2005 ITYPL(IPL)=4 
VPL(IPL)=47 
GO TO 2000 

2006 ITYPL(IPL)=4 
VPL(IPL)=55 
GO TO 2000 

2007 ITYPL(IPL)=1 
VPL(IPL)=0.02D0 
GO TO 2000 

2008 ITYPL(IPL)=1 
VPL(IPL)=0.03D0 
GO TO 2000 

2009 ITYPL(IPL)=1 
VPL(IPL)=0.05D0 
GO TO 2000 

2010 ITYPL(IPL)=1 
VPL(IPL)=0.08D0 
GO TO 2000 

2000 CONTINUE 

C COMPUTE »'PLUME" POINTS AT Y=YN 
DO 1000 IPL=1,KPL 
ITYP=ITYPL(IPL) 

IF(ITYP.EQ.O) GO TO 1000 
GO TO (1,2, 3, 4, 5), ITYP 

1 CONTINUE 

C BREAKDOWN SURFACE PLUME. 

C NOTE THAT DUE TO DIFFERENCE-CENTERING OF 
C Y-COORDINATE IS 0 . 55((YF+YN) , RATHER THAN 
C IN THE PLOTTING CODE. 

BO=VPL(IPL) 

XTEMP(1)=XN(1) 

DO 11 1=2, KN 

XTEMP(I)=0.5D0X(XN(I)+XN(I-D) 

11 CONTINUE 
1=2 

CALL INTERX(1,I,B0,BN,KN,XB0,XTEMP) 

XPL(J,IPL)=XBO 
GO TO 1001 

2 CONTINUE 

C FIND BY INTERPOLATION THE X-COORDINATE WHERE M=MPL 
IF(J.GT.l) GO TO 200 
XPL(J,IPL)=XC 
GO TO 1001 
200 CONTINUE 

MPL=VPL(IPL) 

I=-KN 



GRADIENTS, THE 
YN. IT CAN BE 



ACCURATE 

ADJUSTED 
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CALL INTERX(l,I,MPL,MN,KN,XMO,XN) 

XPL(J,IPL)=XMO 
GO TO 1001 

3 CONTINUE 

C STREAMLINE INTERPOLATION. 

IF(J.GT.l) GO TO 300 
XPL(J,IPL)=VPL(IPL) 

GO TO 1001 

300 CONTINUE 
XSF=XPL(J-1,IPL) 

ISF=2 

ISN=2 

CALL INTERP(0,ISF,KF,XSF,XF,RMSF,RMF,RPSF,RPF) 

CALL RFUNC(RMSF,RPSF,MSF,MUSF,TETASF) 
XSN=XSF+DYXDTAN(PAI2-TETASF) 

ITER=1 

301 ITER=ITER+1 

CALL INTERP(1,ISN,KN,XSN,XN,RMSN,RMN,RPSN,RPN) 

CALL RFUNC(RMSN,RPSN,MSN,MUSN,TETASN) 

TETAAV=0 . 5D0X(TETASF+TETASN) 

XSN=XSF+DY*DTAN(PAI2-TETAAV) 

IF(ITER.LT.ITER0+2) GO TO 301 

XPL(J,IPL)=XSN 

GO TO 1001 

4 CONTINUE 

C CHARACTERISTIC LINE. 

KC=IDINT(VPL(IPL)+l.D-5) 

IF(J.GT.l) GO TO 41 
XPL(J,IPL)=XCHARF(KC) 

GO TO 1001 
41 CONTINUE 

XPL(J,IPL)=XCHARN(KC) 

IFCCSIGNNCKC) .EQ.O.) XPL( J, IPL)=1 . E33 
GO TO 1001 

5 CONTINUE 

C CONSTANT LATERAL (X) OPACITY 
CALL OPACX 
XIC=VPL(IPL) 

DO 51 11=2, KF 
I1=KF-II+1 
12 = 11+1 

XI1=XI(I1, JXI) 

XI2=XI(I2,JXI) 

IF((XIC-XIl)x(XIC-XI2).GT.O.) GO TO 51 
F2=(XI2-XIC)/(XI2-XI1) 

Fl=l .D0-F2 

IF(Fl.LT.O.) CALL FIN(1351) 

IF(F2.LT.O.) CALL FIN(1352) 

XIFC=F25(XF(I1) + F15(XF(I2) 

GO TO 52 

51 CONTINUE 
XIFC=1.D30 

52 CONTINUE 
XPL(J,IPL)=XIFC 
GO TO 1001 

1001 CONTINUE 
IF(J.GT.l) GO TO 1002 
YPL(J)=YC 

GO TO 1000 

1002 CONTINUE 
YPL(J)=YN 

1000 CONTINUE 

SUBTTUUTTTTEnmTJN 

C SUBROUTINE NUMBER 14 

IMPLICIT REAL^8(A-H,L-Z,$) 

REALJ(4 XPL,YPL 

COMMON /PLUME/XPL(1002, 10),YPL(1002) 

COMMON /IPLUME/KPL, ITYPL(IO) 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1 TETAF(101),BF(101), 
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2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), JET1153 

3 TETAN(101),BN(101),XTEMP(101) JET115A 

COMMON/THICKY/XTHd 002), TH( 1002) JET1155 

COMMON /GAMA/G,G1,G2,G3,GA,G5,G6,G7,G8,G9,G10,G11,G12,G13,G1A,G15, JET1156 

1 G16,G17,G18,G19,G20 JET1157 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, JET1158 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, JET1159 

2 TETSYM,TETLIM,DDY,DYMAX JET1160 

COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 JET1161 

COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, JET1162 

1 KF2,IDEL,JDEL, JYXI, JXI,ILEAD,ILEADF,KCLEAD JET1163 

COMMON /ROW/YF,YN,DXF,DXN JET116A 

COMMON /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), JET1165 

1 RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), JET1166 

2 TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), JET1167 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), JET1168 

A MCHARK92) JET1169 

COMMON /ICHARA/KCHARP,KCHARM,KCHARO JET1170 

DIVIDE LINE Y=YN INTO KN-1 INTERVALS. JET1171 

THE X-GRID IS NON-UNIFORMLY DEFINED AS FOLLOWS: JET1172 

(1) (XCHARN(I),YCHARN(D), (XCHARF( I ) , YCHARFC I ) ) , 1 = 1, 2, . . . , KCHARP, JET1173 

DENOTE NEW AND OLD (FORMER) CHARACTERISTIC (C+) POINTS. LET 1=1 JET1174 
AND I=KCHARP CORRESPOND TO THE LEADING AND BOUNDARY JET1175 

CHARACTERISTICS (C+) . JET1176 

(2) THE GRID CONSISTS OF TWO SEGMENTS. THE SO-CALLED FLAT SEGMENT JET1177 

IS BETWEEN X=0 AND X=XLEAD=XCHARN(KCLEAD) . THE SECOND IS THE JET1178 

FAN SEGMENT. IT IS FROM XLEAD TO XBOUND=XCHARN( KCHARP) . JET1179 

(3) THE FAN SEGMENT IS INITIALLY DIVIDED INTO FRACGx(KFO-l) INTERVALSJET1180 



(4) 



(5) 



DEFINED BY THE FAMILY OF C+ CHARACTERISTIC LINES MCHARI(l) TO 
MCHARKKCHARP) . 

THE FLAT SEGMENT IS DIVIDED INTO (l-FRACG)x(KFO-l) EQUAL 
INTERVALS, AS LONG AS THEY ARE NOT SMALLER THAN THE AVERAGE 
FAN INTERVAL. WHEN THEY ARE, THEIR NUMBER IS REDUCED, BUT NOT 
BELOW THREE. 

KCLEAD IS INITIALLY 1. IT IS UPDATED SO THAT THE FLAT SEGMENT 
IS AT LEAST TWICE THE AVERAGE FAN INTERVAL. 

ILEADF=ILEAD 

KCLEAD=0 

DO 1 KC=1, KCHARP 
IF(CSIGNN(KC).LT.O.) GO TO 1 
KCLEAD=KC 

KFAN=KCHARP-KCLEAD 
XL EAD=XCHARN( KCL EAD) 

XBOUND=XCHARN(KCHARP) 

DXl=(XBOUND-XLEAD)/DFLOAT(KFAN) 

IF(XLEAD/DX1 .GT.2.D0) GO TO 11 

1 CONTINUE 
11 CONTINUE 

IF(KCLEAD.EQ. 0) CALL FIN(1401) 

IF(KCLEAD.EQ. KCHARP) CALL FIN(1402) 

ILEAD=IDINT(XLEAD/DXl)+2 
IFCILEAD+KFAN.GT.KFO) I LEAD=KFO-KFAN 
ILEAD1=ILEAD-1 
KN=ILEAD+KFAN 

IF(KN.GT.KFO) CALL FIN(1411) 

DX=XLEAD/DFL0AT(ILEAD1) 

XN(1)=0. 

DO 2 I=1,ILEAD1 
XN(I)=XN(l) + DX5«DFLOAT(I-l) 

2 CONTINUE 

DO 3 I=ILEAD,KN 
XN(I)=XCHARN(KCLEAD+I-ILEAD) 

3 CONTINUE 
RETURN 
END 



ysr£P 



SUB'ROUIINE YSTEP 
C SUBROUTINE NUMBER 15 

IMPLICIT REAL5«8(A-H,L-Z,$) 

REAL3(4 XPL,YPL 

COMMON /PLUME/XPL( 1002, 10 ) ,YPL( 1002) 

COMMON /IPLUME/KPL,ITYPL(10) 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 



JET1181 

JET1182 

JET1183 

JET1184 

JET1185 

JET1186 

JET1187 

JET1188 

JET1189 

JET1190 

JET1191 

JET1192 

JET1193 

JET1194 

JET1195 

JET1196 

JET1197 

JET1198 

JET1199 

JET1200 

JET1201 

JET1202 

JET1203 

JET1204 

JET1205 

JET1206 

JET1207 

JET1208 

JET1209 

JET1210 

JET1211 

JET1212 

JET1213 

JET1214 

JET1215 

JET1216 

JET1217 



TETTZTZ" 

JET1219 

JET1220 

JET1221 

JET1222 

JET1223 

JET1224 
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C 

C 

C 

C 



1 TETAF(101),BF(101), JET1225 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), JET1226 

3 TETAN(101),BN(101),XTEMP(101) JET1227 

COMMON/THICKY/XTH(1002),TH(1002) JET1228 

COMMON /GAMA/G,G1,G2,G3,GA,G5,G6,G7,G8,G9,G10,G11,G12,G13,G1A,G15, JET1229 

1 G16,G17,G18,G19,G20 JET1230 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, JET1231 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, JET1232 

2 TETSYM,TETLIM,DDY,DYMAX JET1233 

COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 JET123A 

COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, JET1235 

1 KF2,IDEL,JDEL, JYXI,JXI,ILEAD,ILEADF,KCLEAD JET1236 

COMMON /ROW/YF,YN,DXF,DXN JET1237 

COMMON /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), JET1238 

1 RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), JET1239 

2 TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), JET1240 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), JET12A1 

A MCHARK92) JET12A2 

COMMON /ICHARA/KCHARP,KCHARM,KCHARO JET12A3 

COMPUTE NEXT Y-STEP. JET12AA 

DYNEXT IS DEFINED AS THE MINIMAL "TRIANGULATION" Y-STEP DYT, 0BTAINEDJET1245 
BY FORWARD INTERSECTION OF C-,C+ CHARACTERISTICS FROM ADJACENT GRID JET12A6 



POINTS XI, X2, 

DYMIN=1.DA0 
DO 1 1=3, KF 
X1=XF(I-1) 

X2=XF(I) 

DX=X2-X1 

TP1=DTAN(TETAF( I-1)-MUF( I-l) ) 

TP2 = DTAN(TETAF(I)+MUF(D) 

F1=-TP2/(TP1-TP2) 

IF(F1.LE.0.) PRINT 555, I , XI , X2, DX, TPl , TP2, FI 
555 F0RMAT(/1X, •I,X1,X2,DX,TP1,TP2,F1=*,I5,6D1A.6/) 
IF(F1.LT.0.) CALL FIN(1501) 

DYT=F1XDXXTP1 

IF(DYT.LE.O.) CALL FIN(1502) 

DYMIN=DMIN1 ( DYMIN, STAB^DYT) 

1 CONTINUE 

DYNEXT=DYMIN 

MOi/e 
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JET12A8 
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JET1255 

JET1256 

JET1257 

JET1258 

JET1259 

JET1260 

JET1261 

JET1262 

JET1263 

JET126A 

JET1265 



SUBROUTINE MOVE 

SUBROUTINE NUMBER 16 

IMPLICIT REALX8(A-H,L-Z,$) 

REAU(A XPL,YPL 

COMMON /PLUME/XPL ( 1002, 10),YPL( 1002) 

COMMON /IPLUME/KPL,ITYPL(10) 

COMMON /VECS/XF(101),RMFa01),RPF(101),MF(101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3 TETAN(101),BN(101),XTEMP(101) 
COMMON/THICKY/XTH(lOO2),TH(10O2) 

COMMON /GAMA/G,G1,G2,G3,GA,G5,G6,G7,G8,G9,G10,G11,G12,G13,G1A, 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI,PAI2, DEG, XC,YC,MEXIT,MFIN,YMAX,DYO,DY, DYNEXT, 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /STAG/RH00,N0,P0,T0,A0,MD0T1 
COMMON /IPAR/JMAX,KFO, ITERO,KF,KN,IM,IP,J, 

1 KF2,IDEL, JDEL, JYXI, JXI,ILEAD,ILEADF,KCLEAD 

COMMON /ROW/YF,YN,DXF,DXN 

COMMON /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), 

1 RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), 

2 TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 

A MCHARK92) 

COMMON /ICHARA/KCHARP, KCHARM, KCHARO 
STORE NEW LINE (N) IN OLD LINE (F). 

KF=KN 

KF2=2XKF 

YF=YN 

DO 1 1=1, KN 
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XFd)=XNd) 


JET1297 


RMFd)=RMNd) 


JET1298 


RPFd)=RPNd) 


JET1299 


MFd)=MNd) 


JET1300 


MUFd)=MUNd) 


JET1301 


TETAFd)=TETANd) 


JET1302 


BFd)=BNd) 


JET1303 


CONTINUE 


JET1304 


DO 2 KC=1,KCHAR0 


JET1305 


IF(CSIGNN(KC).EQ.O.) GO TO 2 


JET1306 


XCHARF(KC)=XCHARN(KC) 


JET1307 


YCHARF( KC) =YCHARN( KC) 


JET1308 


RMCARF(KC)=RMCARN(KC) 


JET1309 


RPCARF(KC)=RPCARN(KC) 


JET1310 


TCHARF(KC)=TCHARN(KC) 


JET1311 


MUCARF( KC) =MUCARN( KC) 


JET1312 


MCHARF(KC)=MCHARN(KC) 


JET1313 


CSIGNF(KC)=CSIGNN(KC) 


JET1314 


CONTINUE 


JET1315 


RETURN 

END 


^ 



SUBROUTINH OHA’CJ? JET1318 

C SUBROUTINE NUMBER 17 JET1319 

IMPLICIT REALX8(A-H,L-Z,$) JET1320 

REALX4 XPL,YPL JET1321 

COMMON /PLUME/XPL(1002,10),YPL(1002) JET1322 

COMMON /IPLUME/KPL,ITYPU10) JET1323 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), JET1324 

1 TETAF(101),BF(101), JET1325 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), JET1326 

3 TETAN(101),BN(101),XTEMP(101) JET1327 

COMMON/THICKY/XTH( 1002) ,TH( 1002) JET1328 

REALX4 YXI,XI,XIPM,XIGRP,XIAPP,XIF JET1329 

COMMON /THICKX/YXI(20),XI(101,20),XIPMa01,20),XIGRP(101,20) JET1330 

1 ,XIAPP(101,20),XIF(101,20) JET1331 

COMMON /GAMA/G,G1,G2,G3,GA,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, JET1332 

1 G16,G17,G18,G19,G20 JET1333 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, JET1334 

1 STAB, DELTA,PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, JET1335 

2 TETSYM,TETLIM,DDY,DYMAX JET1336 

COMMON /STAG/RHOO,NO,PO,TO,AO,MDOT1 JET1337 

COMMON /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), JET1338 

1 RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), JET1339 

2 TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), JET1340 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), JET1341 

4 MCHARK92) JET1342 

COMMON /IPAR/JMAX,KFO, ITERO,KF,KN,IM,IP, J, JET1343 

1 KF2,IDEL, JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD JET1344 

COMMON /ROH/YF,YN,DXF,DXN JET1345 

C COMPUTE X-OPACITY. JET1346 

C BEGIN FROM LIMITING CHARACTERISTIC OF AN ASSUMED P.M. FAN. JET1347 

C XIO — THE THICKNESS BETWEEN THE LIMITING CHARACTERISTIC AND THE JET1348 

C BOUNDARY CHARACTERISTIC OF THE NUMERICAL COMPUTATION. JET1349 

DO 12 1=1, KFO JET1350 

XIFd, JXI) = XF(I) JET1351 

XI(I,JXI)=0. JET1352 

XIPM(I,JXI)=0. JET1353 

XIGRPd, JXI) = 0. JET1354 

XIAPPd,JXI) = 0. JET1355 

12 CONTINUE JET1356 

IF(J.EQ.l) GO TO 1000 JET1357 

PSILIM=TETLIM JET1358 

XLIM=XC+(YF-YC)/DTAN(PSILIM) JET1359 

XBOUND=XF(KF) JET1360 

KPM=10 JET1361 

DX=(XLIM-XBOUND)/DFLOAT(KPM) JET1362 

SUM=0. JET1363 

DO 1 1=1, KPM . JET1364 

X1=XB0UND+DFL0AT(I-1)XDX JET1365 

X2=X1+DX JET1366 

PSl=PAI2-DATAN((Xl-XC)/(YF-YO) JET1367 

PS2=PAI2-DATAN( (X2-XC)/(YF-YC) ) JET1368 
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16 

15 

14 



21 

2 

C 
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Q1=(PS1-PSILIM)/G5 
Q2=(PS2-PSILIM)/G5 
IF(I.EQ.KPM) Q2=1.D-10 
IF(Q2.LT.O.) CALL FIN(1701) 

F1=G11x(DSIN(Q 1))XX(2.DO/(G-1 .DO)) 
F2=GllX(DSIN(Q2))XX(2.DO/(G-l.DO)) 

SUM=SUM+DX*(Fl+F2)/2.DO 

CONTINUE 

XI0=SUMX(N0XSIGMA) 

RE-EVALUATE XIO FOR A RING-JET. 

IF(DELTA.EQ.O.) GO TO 14 
M=MFIN 

CALL MFUNC(M,F,ETA,TETA) 

PSI=TETA+DARSIN(l.DO/M) 

GOREM=l .DO+G1XMXX2 
G0R=MXX2-1 .DO 
CALL HINTER(M,HM) 

DELT0B=0.5D0XDSQRT(G0R)X(l ,DO/(MEXITxETA)+DSIN(TETA)/M)/DSIN(PSI) 
1 +G15XHM/2.DO 

EVER=SIGMAXN0XYC/(MXDSIN(TETA)*DSIN(PSI)XG0REMXXG6) 
GGG=2.DO-DELTOBX(G+l .DO)/2.DO 
IFCDABS(GGG) .GT.l.D-10) GO TO 15 
PRINT 16, DELTOB,G,GGG 

FORMATC/IX, 'FROM OPACX. GGG NEARLY ZERO. EXPRESSION FOR XIO IS', 

1 IX, 'SINGULAR. DELT0B,G,GGG=',3D12.4/) 

CALL FIN(1715) 

CONTINUE 

EVER=EVER/GGG 

XIO = EVER3<( (YF/YOXXGGG-l . D0)/(YF/YC) 

CONTINUE 

XI(KF,JXI)=XIO 

XIPM(KF,JXI)=XIO 

XIGRP(KF,JXI)=XIO 

KF1=KF-1 

DO 2 11=1, KFl 

I=KF-II+1 

X1=XF(I) 

X2=XF(I-1) 

DX=X1-X2 

F1 = 1.D0/(1.D0+G1XMF(I )5<X2)XXG6 

F2=l.D0/(l .D0+G1XMF(I-1)XX2)XXG6 
DTNUM=(N0XSIGMA)XDXX(F1+F2)/2.D0 
XKI-l, JXI)=XI(I, JXD + DTNUM 
XIPM(I-1, JXI)=1 .D24 
XIGRP(I-1,JXI)=1.D24 
PS1=PAI2-DATAN( (Xl-XC)/(YF-YO) 

PS2=PAI2-DATAN( (X2-XC)/(YF-YC) ) 

IF(PS2.GT.PSI1) GO TO 2 
Q1=(PS1-PSILIM)/G5 
Q2=(PS2-PSILIM)/G5 
IF(Q1.LT.0.) CALL FIN(1711) 

F1=G11X(DSIN(Q1))XX(2.D0/(G-1 .DO)) 
F2=G11X(DSIN(Q2))XX(2.D0/(G-1.D0)) 

DTPM=(N05(SIGMA)XDXX(F1 + F2)/2.D0 
XIPM(I-1, JXI)=XIPM(I, JXD + DTPM 
DISTl=DSQRT((Xl-XC)X5f2+(YF-YC)xx2) 
DIST2=DSQRT((X2-XC)X3«2+(YF-YC)5<5<2) 

KC1=KCLEAD+I-ILEAD 
KC2=KC1— 1 

IF(KC2.LT.KCLEAD) GO TO 21 
M1=MCHARI(KC1) 

M2=MCHARI(KC2) 

CALL MATCHCI , Ml , MGl ,M0BI1 ,MABI 1 ) 

CALL MATCH(I-1,M2,MG2,M0BI2,MABI2) 

Fl=l .DO/(l .D0+G1XMG1X3(2)XXG6 
F2=l .D0/(1 .D0+G1XMG2XX2)3(XG6 
DTGRP=(N0XSIGMA)3<DX3<(F1 + F2)/2.D0 
XIGRP(I-1,JXI)=XIGRP(I, JXD+DTGRP 
CONTINUE 
CONTINUE 

APPROXIMATE THICKNESS XIAPP( I , JXI ) . BASED ON CLOSED-FORM INTEGRATION 
DO 3 1=1, KF 
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XIAPPd, JXI) = 1,D2A 
KC=KCLEAD+(I-ILEAD) 

IFCDELTA.EQ.O.) GO TO 3 
IF(KC.LT.KCLEAD) GO TO 3 
IF(XF(I) .LT.XCHARF(D) GO TO 3 
M=MCHARI(KC) 

CALL MFUNC(M,F,ETA,TETA) 

PSI=TETA+DARSIN( 1 .DO/M) 

GOREM=1.DO+G1XM3(?€2 
G0R=M*X2-1 .DO 
CALL HINTER(M,HM) 

DELTOB=0.5DOXDSQRT(GOR)3<(1 .D0/(MEXITXETA)+DSIN(TETA)/M)/DSIN(PSI) 
1 +G15XHM/2.DO 

EVER=SIGMA5«NOXYC/(M?(DSIN(TETA)5(DSIN(PSI)5(GOREM5(XG6) 

GGG=2 . DO-DELTOBx( G+1 . DO )/2 . DO 
IF(DABS(GGG) .GT.l.D-10) GO TO 25 
PRINT 26, I,KC,M,DELTOB,G,GGG 

26 F0RMAT(/1X, 'FROM OPACX. GGG NEARLY ZERO. EXPRESSION FOR XIO IS', 

1 IX, 'SINGULAR. I,KC,M=',I5,D12.A/ 

2 IX, 'DELT0B,G,GGG=',3D12.4/) 

CALL FIN(1725) 

25 CONTINUE 

EVER=EVER/GGG 

XIAPPd, JXI) = EVER*((YF/YC)X3(GGG-1.D0)/(YF/YC) 

3 CONTINUE 
1000 CONTINUE 
RETURN 
END 
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SUBROUTINETOTnJC 

SUBROUTINE NUMBER 18 

IMPLICIT REALX8(A-H,L-Z,$) 

REAL5(4 XPL,YPL 

/PLUME/XPL (1002, 10),YPLd 002) 

/IPLUME/KPL,ITYPLdO) 

/VECS/XFd01),RMFd01),RPFd01),MFd01),MUFd01), 
TETAFd01),BFd01), 

XN(101),RMNd01),RPNd01),MN(101),MUNd01), 
TETANd01),BN(101),XTEMP(101) 
COMMON/THICKY/XTHd002),THd002) 

REAL5<4 YXI,XI,XIPM,XIGRP,XIAPP,XIF 

COMMON /THICKX/YXI(20),XId01,20) ,XIPM(1 01,20 ),XIGRP( 101, 20) 

L ,XIAPP(101,20),XIF(101,20) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15. 
L G16,G17,G18,G19,G20 

COMMON /PAR/PAI , PAI2 , DEG , XC, YC, MEXIT, MFIN, YMAX, DYO , DY, DYNEXT , 

STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 
TETSYM, TETLIM, DDY, DYMAX 
/STAG/RHOO,NO,PO,TO,AO,MDOT1 
/IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

KF2,IDEL,JDEL, JYXI,JXI,ILEAD,ILEADF,KCLEAD 
/ROW/YF,YN,DXF,DXN 

/CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), 
RMCARFC92) , RPCARFC92) , RMCARN(92) , RPCARNC92) , 
TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), 
CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), 
MCHARK92) 

/ICHARA/KCHARP,KCHARM,KCHARO 

VARIABLES OF GRID POINTS IN THE FAN SEGMENT FROM THE 
SEMI-INVERSE INTEGRATION (IN SUBR. SEMINV) . NOTE THAT GRID POINTS 
XNd) WERE ALREADY DETERMINED IN SUBR. GRIDN. 

DO 1 I=ILEAD,KN 
KC=KCLEAD+I-ILEAD 
IF(KC.GT.KCHARP) CALL FINdSOl) 

RMN(I)=RMCARN(KC) 

RPN(I)=RPCARN(KC) 

MN(I)=MCHARN(KC) 

MUN(I)=MUCARN(KC) 

TETAN(I)=TCHARN(KC) 

CONTINUE 

- hluFut^c 



c 

c 

c 



COMMON 

COMMON 

1 

COMMON 

COMMON 

1 
2 

3 

4 

COMMON 
LOAD FLOW 



DOUBLE PRECISION FUNCTION NUFUNC(M) 
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C SUBROUTINE NUMBER 19 JET1513 

IMPLICIT REALX8(A-H,L-Z,$) JET1514 

REALXA XPL,YPL JET1515 

COMMON /PLUME/XPL(1002,10),YPL(1002) JET1516 

COMMON /IPLUME/KPL,ITYPL(10) JET1517 

COMMON /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), JET1518 

1 TETAF(101),BF(101), JET1519 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), JET1520 

3 TETAN(101),BN(101),XTEMP(101) JET1521 

COMMON/THICKY/XTH(1002),TH(1002) JET1522 

REALJ«4 YXI,XI,XIPM,XIGRP,XIAPP,XIF JET1523 

COMMON /THICKX/YXI(20),XI(101,20),XIPM(101,20),XIGRP(101,20) JET1524 

1 ,XIAPP(101,20),XIF(101,20) JET152S 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,JET1526 

1 G16,G17,G18,G19,G20 JET1527 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, JET1528 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, JET1529 

2 TETSYM,TETLIM,DDY,DYMAX JET1530 

COMMON /STAG/RHOO,NO,PO,TO, AO,MDOT1 JET1531 

COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, JET1532 

1 KF2,IDEL,JDEL, JYXI,JXI,ILEAD,ILEADF,KCLEAD JET1533 

COMMON /ROW/YF,YN,DXF,DXN JET1534 

COMMON /CHARAC/XCHARF(92),YCHARF(92),XCHARN(92),YCHARN(92), JET1535 

1 RMCARF(92),RPCARF(92),RMCARN(92),RPCARN(92), JET1536 

2 TCHARF(92),TCHARN(92),MUCARF(92),MUCARN(92), JET1537 

3 CSIGNN(92),CSIGNF(92),MCHARN(92),MCHARF(92), JET1538 

4 MCHARK92) JET1539 

COMMON /ICHARA/KCHARP,KCHARM,KCHARO JET1540 

C COMPUTE NU AS FUNCTION OF MACH NUMBER M. NOTE THAT THE P.M. JET1541 

C DEFINITION OF NU HAS BEEN MODIFIED BY ADDING A CONSTANT. THE USUAL JET1542 

C CHOICE OF THE CONSTANT IS SUCH THAT NU=0 FOR INFINITE M. JET1543 

Q=1.D0/DSQRT(MXX2-1.D0) JET1544 

NUFUNC=NU0-(G5XDATAN(G5XQ)-DATAN(Q)) JET1545 

RETURN U.kjt'^CT' JET1546 

END HMOC. I JET1547 

SUBROUTINE HM5EI JET154'8’ 

C SUBROUTINE NUMBER 20 JET1549 

IMPLICIT REAU«8(A-H,L-Z,$) JET1550 

REALMS KAPAOB JET1551 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, JET1552 

1 G16,G17,G18,G19,G20 JET1553 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, JET1554 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, JET1555 

2 TETSYM,TETLIM,DDY,DYMAX JET1556 

COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, JET1557 

1 KF2,IDEL,JDEL,JYXI, JXI,ILEAD,ILEADF,KCLEAD JET1558 

COMMON /GRP/DMINV,MHINV(101),HMV(101) JET1559 

COMMON /IGRP/KHM JET1560 

C A ROUTINE FOR THE C+ DERIVATIVE DUE TO RING SYMMETRY (GRP). JET1561 

KHM=51 JET1562 

IF(KHM.GT.lOl) CALL FIN(2001) JET1563 

MINV0=1 .DO/MEXIT JET1564 

DMINV=MINVO/DFLOAT(KHM-1) JET1565 

M=MEXIT JET1566 

SUM=0. JET1567 

KHM1=KHM-1 JET1568 

DO 1 I=1,KHM1 JET1569 

MF=M JET1570 

MHINV(I)=MINVO-DFLOAT(I-1)5(DMINV JET1571 

M=1.D0/MHINV(I) JET1572 

DM=M-MF JET1573 

M1=M-DM JET1574 

M2=M-DM/2.D0 JET1575 

M3=M JET1576 

CALL MFUNC(M1,F1,ETA1,TETA1) JET1577 

CALL MFUNC(M2,F2,ETA2,TETA2) JET1578 

CALL MFUNC(M3,F3,ETA3,TETA3) JET1579 

SUM=SUM+DMX(Fl+4.D0XF2+F3)/6 .DO JET1580 

ETA=ETA3 JET1581 

TETA=TETA3 JET1582 

PSI=TETA+DARSIN(1.D0/M) JET1583 

NORM=( (3 . D0-G)/4 . D0)X(MXX2-1 . D0)x*0 .75D0/ JET1584 
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1 (DSIN(PSI)*(1.D0+G1*MK*2)XXG1A) 

HM=SUM*NORM 

GOREM=1.DO+G1*M**2 



11 

12 

1 



GOR=MXX2-1.DO 

DELT0B=0.5D0XDSQRT(GOR)X(l.D0/(MEXIT*ETA)+DSIN(TETA)/M)/DSIN(PSI) 

1 +((G+1.D0)/(2.D0X(3.D0-G)))XHM 

EPSIOB=DELTOB/DSQRT(GOR)-DSIN(TETA)/(MXDSIN(PSD) 

KAPA0B=1.D0 

IF(DABS(PAI2-TETA) .GT.l .D-6) 

1KAPA0B=DTAN(TETA)*EPSI0B 
LAMDOB=EPSIOB-DELTOBXGOREM/(GOR*DSQRT(GOR)) 

PRINT 11,I,M,HM,TETAXDEG,PSI*DEG 

F0RMAT(/1X, • I, M,HM,TETA,PSI=», 15,5012.4) 

PRINT 12,DELT0B,EPSI0BXDEG,KAPA0BXDEG,LAMD0B*DEG 
FORMAT( IX, 'DELT0B,EPSI0B,KAPA0B,LAMD0B=',5X,5D12.4) 

CONTINUE 

MHINV(KHM)=0. 

HMV(KHM)=1.D0 

iS""" MFUNC 

SUBROUTINE' HFDNC(M,F,tlA, I tl A) 

SUBROUTINE NUMBER 21 

IMPLICIT REALX8(A-H,L-Z,$) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI, PAI2, DEG,XC, YC, MEXIT,MFIN,YMAX, DYO, DY, DYNEXT, 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 NUPT1,TETLIM 

COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

1 KF2,IDEL, JOEL, JYXI, JXI,ILEAD,ILEADF,KCLEAD 

COMMON /GRP/DMINV,MHINV(101),HMV(101) 



NU=NUFUNC(M) 

TETA=NUFUNC(MEXIT)+PAI2-NU 
G0REM=1 .D0+G1XMXX2 



G0R=MXX2-1 .DO 

F=(M**2)X(GOREMX*G13)XDSIN(TETA)/GORXX1.25DO 
G0REM1=1 . D0+G1*MEXIT**2 
G0R1=MEXITXX2-1 .DO 

ETA=((GOREM/GOREM1)**G14)*((GOR1/GOR)**0 .25D0) 

RETURN 

END 

SUBROUTINE HINTER(M,H) 

C SUBROUTINE NUMBER 22 

IMPLICIT REALX8(A-H,L-Z,$) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI, PAI2, DEG, XC,YC,MEXIT,MFIN,YMAX, DYO, DY, DYNEXT, 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /IPAR/ JMAX, KFO , ITERO , KF, KN, IM, IP, J , 

1 KF2,IDEL, JDEL, JYXI, JXI,ILEAD,ILEADF,KCLEAD 

COMMON /GRP/DMINV,MHINV(101),HMV(101) 

COMMON /IGRP/KHM 
C COMPUTE H(M) BY INTERPOLATION 
MINV=1.D0/M 

I=KHM-IDINT(MINV/DMINV-l.D-9)-l 
IF(I.GE.l.AND.I.LT.KHM) GO TO 1 
PRINT 11,I,KHM,M,MEXIT 

11 F0RMAT(/1X, 'I,KHM,M,MEXIT=',2I5,2D14.6/) 

CALL FIN(2201) 

1 CONTINUE 

F1=(MINV-MHINV(I+1))/DMINV 
F2=l .DO-Fl 



IF(Fl.LT.-l.D-9) CALL FIN(2210) 
IF(F2.LT.-1 .D-9) CALL FIN(22I1) 
H=F15(HMV(I)+F2*HMV(I+1) 

RETURN 

END 

SUBROUTINE MATCHC I , MOB, MAB, MOBI ,MABI ) 
C SUBROUTINE NUMBER 23 
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IMPLICIT REAL5«8(A-H,L-Z,$) 

COMMON /VECS/XF(101),RMFa01),RPFa01),MF(101),MUF(101), 

1 TETAF(101),BF(101), 

2 XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3 TETAN(101),BN(101),XTEMP(101) 

COMMON /ROW/YF,YN,DXF,DXN 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,GI3,GIA,GI5 
1 G16,G17,G18,G19,G20 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DYO,DY,DYNEXT, 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2 TETSYM,TETLIM,DDY,DYMAX 
COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

1 KF2,IDEL, JDEL, JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON /GRP/DMINV,MHINVa01),HMV(101) 

COMMON /IGRP/KHM 

C COMPUTE H(M) AND THE ALFA-DERIVATIVES 
M=MOB 

CALL MFUNC(M,F,ETA,TETA) 

PSI=TETA+DARSIN(1 .DO/MI 
CALL HINTER(M,HM) 

GOREM=1.DO+G15<MXX2 

G0R=MXX2-1.D0 

DELT0B=0 .5DOXDSQRT(GOR)X(l.DO/(MEXIT*ETA)+DSIN(TETA)/M)/DSIN(PSn 
1 +G155<HM/2.D0 

F0B=(G7XG0REM)xxG2/M 
FAB = FOBX(YF/YC)X)(DELTOB 
CALL AREAF(FAB,MAB) 

C COMPUTE MABI FROM THE INVERSE PROBLEM SOLUTION 
COTAV=(XF(I)-XC)/(YF-YC) 

PSI0*PAI2-DATAN(C0TAV) 

EVY=YFXDL0G(YF/YC)/(YF-YC)-I .DO 

PSIN=PSIO 

DO 1 ITER=1,50 

PSI=PSIN 

M=DSQRTa.DO+GA/DTAN((PSI-TETLIM)/G5)XX2) 

M=DMAXI(M,MEXIT) 

CALL HINTER(M,HM) 

CALL MFUNC(M,F,ETA,TETA) 

G0REM=I.D0+GI^MX5(2 
G0R=MXX2-1 .DO 

DELT0B=0.5D0^DSQRT(GOR)x(l.D0/(MEXIT5fETA)+DSIN(TETA)/M)/DSIN(PSI) 
1 +G155«HM/2.D0 

EPSIOB = DELTOB/DSQRT(GOR)-DSIN(TETA)/(M5(DSIN(PSI) ) 

LAMDOB = EPSIOB-DELTOB5«GOREM/(GOR5«DSQRT(GOR)) 

COTN=COTAV+LAMDOB?«EVY/DSIN(PSI)5<X2 

PSIN=PAI2-DATAN(C0TN) 

DPSI=PSIN-PSI 

IF(DABS(DPSI) .LT.l.D-9) GO TO 11 

I CONTINUE 

PRINT 12,I,ITER,PSI,PSIN,DPSI,M,XF(I),YF,XC,YC 
12 F0RMAT(/1X, ' I , ITER, PSI , PSIN, DPSI ,M, XF( I) , YF, XC, YC= '// 

1 1X,2I4,8D11.3/) 

CALL FIN(2301) 

II CONTINUE 

C USING MOBI=M AS COMPUTED FROM THE INVERSE PROBLEM, FIND MABI. 

MOBI=M 

M=M0BI 

CALL MFUNC(M,F,ETA,TETA) 

PSI=TETA+DARSIN(l.DO/M) 

CALL HINTER(M,HM) 

GOREM=l .D0+G15«MX5(2 



C 



GOR=M^?«2-1.DO 

DELT0B = 0 .5D05«DSQRT(GOR)X(1.DO/(MEXIT5«ETA) + DSIN(TETA)/M)/DSIN( 



1 +G15JCHM/2.DO 
F0B = (G75«G0REM)XXG2/M 
FAB = F0B5<(YF/YC)X5«DELT0B 
CALL AREAF( FAB, MABI) 

RETURN 

END 

rUBROUTlNE AREAF(F,M) 

SUBROUTINE NUMBER 24 

IMPLICIT REALX8(A-H,L-Z,$) 
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COMMON /GAMA/G,G1,G2,G3,GA,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, JET1729 

1 G16,G17,G18,G19,G20 JET1730 

COMMON /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, JET1731 

1 STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, JET1732 

2 TETSYM,TETLIM,DDY,DYMAX JET1733 

COMMON /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, JET1734 

1 KF2,IDEL, JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD JET1735 

COMMON /GRP/DMINV,MHINV(101),HMV(101) JET1736 

COMMON /IGRP/KHM JET1737 

C COMPUTE MACH NUMBER M FROM AREA RATIO FUNCTION F JET1738 

C F=((2/(G+I))xa+(G-1)XM*X2))XX((G+1)/(2X(G-1)))/M JET1739 

C INITIAL GUESS IS MIN JET1740 

E1=(FxMEXIT)XX(1.D0/G2)/G7 JET1741 

E2=(E1-1.D0)/G1 JET1742 

E3=DMAX1(E2,MEXITXX2) JET1743 

MIN=DSQRT(E3) JET1744 

EMN=MIN JET1745 

DO 1 1=1,100 JET1746 

EMO=EMN JET1747 

G0REM=1.D0+G1*EM03(X2 JET1748 

G0R=EM0XX2-1 .DO JET1749 

F0=(G7XG0REM)XXG2/EM0 JET1750 

DF=FO-F JET1751 

C PRINT 123,I,EM0,EMN,F0,F,DF,G0R,G0REM JET1752 

C123 FORMATCIX, *I,EM0,EMN,F0,F,DF,G0R,G0REM=M5,7D12.4) JET1753 

DFDM=F0XG0R/(EM0XG0REM) JET1754 

DMN=DF/DFDM JET1755 

EMN=EMO-DMN JET1756 

EPSEM=DABS(DMN/EMN) JET1757 

IFCEPSEM.lt. l.D-10) GO TO 11 JET1758 

I CONTINUE JET1759 

CALL FINC2401) JET1760 

II CONTINUE JET1761 

M=EMN JET1762 

RETURN JET1763 

END JET1764 
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